Localization, synchronization and navigation using passive sensor networks

ABSTRACT

A method for sensor operation includes deploying a network of sensors ( 22 ), which have respective clocks ( 36 ) that are not mutually synchronized. At least a group of the sensors receive respective signals emitted from each of a plurality of sources ( 24, 26 ), and record respective times of arrival of the signals at the sensors according to the respective clocks. Location information is provided, including respective sensor locations of the sensors. The respective clocks are synchronized based on the recorded times of arrival and on the location information. In the process the sources may be localized, or if the sources are far away, then their directions may be resolved. Sensor positions may also be resolved in the process.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of U.S. Provisional PatentApplication 61/617,186, filed Mar. 29, 2012, and of U.S. ProvisionalPatent Application 61/755,485, filed Jan. 23, 2013, both of which areincorporated herein by reference.

FIELD OF THE INVENTION

The present invention relates generally to sensor networks, andparticularly to methods and systems for synchronization, localizationand navigation using sensor networks.

BACKGROUND

Source localization systems based on Time of Arrival (TOA) are used tolocate emitters or reflectors using sensors distributed in an area ofinterest. The area can be as large as the earth (as in Global SatelliteNavigation Systems, GNSS, including the well-known GPS) to very small,such as a room in an office building. Each sensor typically receivessignals from emitters in its range, records the times of arrival of thesignals, and then reports the results to a computer (at a singlelocation or distributed), which constructs an estimation of the emitterlocations. In many systems, the signals are processed to produceDifferential Time of Arrival (DTOA) measurements, based on thedifferences between times of arrival of a signal at different locations.For convenience in the description that follows and in the claims, theterm “TOA” is used generically to refer to all forms of measurement ofsignal times of arrival, including DTOA, unless specifically indicatedotherwise. The term “emitter” is used generically to refer to all formsof point energy sources, including bodies that reflect energy.

Sensor synchronization is a key factor in the ability of the computer tofind the emitter locations. Synchronization may be maintained byequipping each sensor with a very accurate clock (such as an atomicclock) or a GPS receiver, but these solutions are costly, pose hardconstraints on system design, and require extra electrical power. Insome systems, synchronization is maintained by two-way messaging ortwo-way TOA ranging between the sensors, but these approaches alsoincrease the cost, complexity and power consumption of the sensors.

Critical infrastructures such as wireless communication networks, stocktrading systems, smart power grids, airport landing guiding systems, anddigital broadcast networks, are heavily dependent on GNSS signals fortime and positioning. Since the GNSS signal is weak, however, even thesmallest jammer can do significant damage to infrastructure, and manyjamming events are detected daily around the world. Furthermore, if theoperator of the GNSS constellation chooses to shut down its services (orperhaps limit them to certain authorized military users, for example),GNSS-dependent infrastructures will fail.

SUMMARY

Embodiments of the present invention that are described hereinbelowprovide methods and systems that use measurements made by a network ofunsynchronized sensors to perform accurate synchronization of the sensorclocks, as well as localization and navigation using this sort ofsynchronization.

There is therefore provided, in accordance with an embodiment of thepresent invention, a method for sensor operation, which includesdeploying a network of sensors, the sensors having respective clocksthat are not mutually synchronized. At least a group of the sensorsreceive respective signals emitted from each of a plurality of sources,and record respective times of arrival of the signals at the sensorsaccording to the respective clocks. Location information is provided,including respective sensor locations of the sensors. The respectiveclocks are synchronized based on the recorded times of arrival and onthe location information.

In a disclosed embodiment, synchronizing the respective clocks includesestimating offsets and skews between the respective clocks.

The method may also include computing the source locations based on thesensor locations and the recorded times of arrival. In a disclosedembodiment, computing the source locations includes applying anestimator to a set of equations relating the recorded times of arrivaland the source and sensor locations. Applying the maximum likelihoodestimator may include applying an iterative optimization process to theset of the equations, wherein the optimization process derives a set oflinear constraints on offsets and skews of the respective clocks basedon the received signals. Additionally or alternatively, applying theiterative optimization process includes performing a convex optimizationusing maximum volume inscribed ellipsoid centering.

In one embodiment, the method includes detecting a fault in the networkof the sensors based on the location information and the recorded timesof arrival.

Typically, receiving the respective signals includes receiving radiosignals, wherein the sources may include at least one satellite source,or possibly three or more satellite sources. Alternatively, theplurality of the sources may include the at least one satellite sourceand at least one terrestrial source. Synchronizing the respective clocksmay then include using the times of arrival of the respective signalsemitted from only a single satellite source and a single terrestrialsource in order to synchronize the respective clocks of the sensors.

In some embodiments, the method includes finding a location of the atleast one terrestrial source based on the recorded times of arrivaland/or finding a direction of the at least one satellite source based onthe recorded times of arrival.

Typically, the at least one satellite source is not a Global SatelliteNavigation Systems (GNSS) satellite. Synchronizing the respective clocksmay include detecting an attempt to spoof a satellite source, anddiscarding the signals received from the spoofed satellite source.

In a disclosed embodiment, synchronizing the respective clocks includesapplying an orthogonal decomposition to a measurement space of thesatellites.

Additionally or alternatively, receiving the respective signalsincludes, upon detecting a loss of signal from one of the sources,selecting a new source, and recording the respective times of arrival ofthe signals from the new source.

In one embodiment, the method includes navigating based on the recordedtimes of arrival and the location information.

In another embodiment, the sensors are associated with respectivecomputers, and the method includes synchronizing operation of thecomputers based on synchronization of the respective clocks.

In an alternative embodiment, receiving the respective signals includesreceiving acoustic signals.

There is also provided, in accordance with an embodiment of the presentinvention, a method for transmitter operation, which includes deployinga network of sources, having respective clocks that are not mutuallysynchronized. At least a group of the sources transmit, at timesdetermined according to the respective clocks, respective signals, whichare received by a plurality of sensors, and record respective times ofarrival of the signals at the sensors. Location information is provided,including respective source locations of the sources. The respectiveclocks are synchronized based on the recorded times of arrival and onthe location information. The sensors may be located accordingly.

In one embodiment, the network of sources includes base stations in acellular communications network, and the sensors include user equipmentin the cellular communications network, and synchronizing the respectiveclocks includes synchronizing operation of the base stations based onthe signals received by the user equipment.

There is additionally provided, in accordance with an embodiment of thepresent invention, a method for sensor operation, which includesdeploying a network of sensors, the sensors having respective clocksthat are not mutually synchronized. At least a group of the sensorsreceive respective signals emitted from each of a plurality of sources,and record respective times of arrival of the signals at the sensorsaccording to the respective clocks. Location information is provided,including respective sensor locations of the sensors. The sourcelocations are computed based on the recorded times of arrival and on thelocation information.

There is further provided, a network system, including a network ofsensors, which include respective clocks that are not mutuallysynchronized, and which are configured to receive, in at least a groupof the sensors, respective signals emitted from each of a plurality ofsources, and to record respective times of arrival of the signals at thesensors according to the respective clocks. A processor is configured toprocess the recorded times of arrival, using location informationincluding respective sensor locations of the sensors, so as tosynchronize the respective clocks.

There is moreover provided, in accordance with an embodiment of thepresent invention, processing apparatus, including a communicationsinterface, which is configured to receive, from at least a group ofsensors in a network of the sensors, which have respective clocks thatare not mutually synchronized, times of arrival recorded by the sensorsaccording to the respective clocks of signals emitted respectively fromeach of a plurality of sources. A processor is configured to process therecorded times of arrival, using location information includingrespective sensor locations of the sensors, so as to synchronize therespective clocks.

There is furthermore provided, in accordance with an embodiment of thepresent invention, a computer software product, including anon-transient computer-readable medium in which program instructions arestored, which instructions, when read by a computer, cause the computerto receive, from at least a group of sensors in a network of thesensors, which have respective clocks that are not mutuallysynchronized, times of arrival recorded by the sensors according to therespective clocks of signals emitted respectively from each of aplurality of sources, and to process the recorded times of arrival,using location information including respective sensor locations of thesensors, so as to synchronize the respective clocks.

Other embodiments provide systems, apparatus and software products thatoperate in accordance with the other methods laid out above.

The present invention will be more fully understood from the followingdetailed description of the embodiments thereof, taken together with thedrawings in which:

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram that schematically illustrates a sensornetwork, in accordance with an embodiment of the present invention;

FIG. 2 is a block diagram that schematically illustrates components of asensor in the network of FIG. 1;

FIG. 3 is a block diagram that schematically illustrates components of asynchronization controller in the network of FIG. 1;

FIG. 4 is a flow chart that schematically illustrates a method for clocksynchronization and localization in a sensor network, in accordance withan embodiment of the present invention;

FIG. 5 is a schematic, pictorial illustration of a navigation systembased on a sensor network, in accordance with an embodiment of thepresent invention;

FIG. 6 is a schematic, pictorial illustration of a computer system withsensor-based synchronization, in accordance with an embodiment of thepresent invention; and

FIG. 7 is a schematic, pictorial illustration of a communication systemwith sensor-based synchronization, in accordance with an embodiment ofthe present invention.

DETAILED DESCRIPTION OF EMBODIMENTS Overview

Embodiments of the present invention that are described herein deriveaccurate synchronization from networks of sensors, without requiringthat the respective clocks of the sensors themselves be mutuallysynchronized. In the disclosed embodiments, a group of sensors receivesrespective signals emitted from multiple sources, which may beterrestrial, space-based (i.e., satellites), or a combination of thetwo. The sensors record the times of arrival of the signals at thesensor locations according to their respective, unsynchronized clocks.

A processor (which may be centralized or associated with one or more ofthe sensors) receives and uses the recorded times of arrival, togetherwith known location information, to synchronize the sensor clocks.Typically, this known location information comprises the sensorlocations, while the source locations are not known a priori, althoughsome of the source locations may be known. “Synchronizing the clocks,”in the context of the present description and in the claims, typicallymeans calculating the offsets and skews between the variousunsynchronized sensor clocks, whether or not this synchronization isactually applied by the sensors themselves or only as a part of thecomputations performed by the processor. In some cases, some of thesensor clocks may be synchronized in advance, while others are not; andthe term “synchronizing the clocks,” as used in the present descriptionand in the claims includes cases in which only some of the clocks needbe synchronized. As a by-product of the synchronization computation, theprocessor may also compute the hitherto-unknown source or sensorlocations.

In alternative embodiments, the same methods of synchronization andlocalization may be applied, mutatis mutandis, to a network ofunsynchronized emitters. Such techniques may be applied, for example, insynchronizing a cellular communication network, as explained in detailhereinbelow. For simplicity and clarity, the methods described belowwill refer in detail to the case of sensor synchronization, but themodification of these methods for emitter synchronization will bestraightforward for those skilled in the art and is considered to bewithin the scope of the present invention.

The processor thus synchronizes and calibrates the time offsets andskews of the sensors (or, in alternative embodiments, the emitters) inthe system passively, based only on the information that the sensorsreceive by wireless reception from source transmissions. As a result,the sensors can use inaccurate internal clocks and need only to becapable of measuring the time of arrival of the emitter signals relativeto their own clocks and to relay this information to the processor. Thesensors can thus be made small and inexpensive, while enjoying very lowpower consumption. The sources are not generally required to cooperatewith the system in any way, thus allowing the use of sources ofconvenience—such as existing satellite or terrestrial sources, or verysimple dedicated signal emitters—which operate in parallel without anyinterconnection.

Only limited location information need be known in advance for operationof the synchronization methods that are described herein. For example,the known locations (of the sensors or the sources) may be given inabsolute terms or only relative to one another. Furthermore, if thelocations of a certain number of the sensors (typically at least three)are already known, the remaining sensors can be located if they areallowed to transmit until their locations are estimated to the requiredaccuracy. Alternatively, the sensors of unknown position can beconsidered navigators, and their positions and clock offsets can beresolved without necessarily ranging transmission on their part.

In some of the embodiments that are described hereinbelow, the processorcarries out a process of estimation and optimization, such as maximumlikelihood estimation, in order to compute both the offsets (and theskews if any) between the respective clocks and the unknown sensorlocations. The possible values of the clock offsets and skews can belinearly constrained by the time of arrival measurements of the sourcesignals. The clock offsets and skews and the unknown locations areestimated by applying an iterative optimization process to a set ofequations and constraints relating the recorded times of arrival and thesource and sensor locations. For terrestrial sources in particular,convex optimization using maximum volume inscribed ellipsoid centeringhas been found to give a good starting point for this estimationprocess. When satellite sources are present, the convex optimization canbe replaced by other means, as described in detail hereinbelow.

Embodiments of the present invention are capable of operating withdifferent numbers and configurations of terrestrial and/or satellitesignal sources. For example, synchronization is possible in all of thefollowing combinations:

-   -   Satellite sources only (at least three sources);    -   Terrestrial sources (at least two, depending on the number of        sensors);    -   Combination of satellite and terrestrial sources (at least one        of each, and only one of each is needed if at least six or seven        sensors are used).        Synchronization can be improved by adding sources, as well as        sensors. The actual positions of the satellite sources in the        sky have only a second-order effect on synchronization and        localization accuracy. In every embodiment in which a sensor        network is synchronized, this synchronization can be used to        localize terrestrial sources, as well as to find the direction        of a satellite or navigate a sensor.

Precision of synchronization in sensor networks according to embodimentsof the present invention is limited only by the available processingcapability and the sensor position accuracy. The precision can beenhanced by increasing the bandwidth, signal/noise ratio, observationtime and number of sources processed. The sensor positions may bedetermined by any suitable means known in the art, such as surveying,localization (if the sensors also transmit signals), or navigation bythe methods described herein, assuming that the positions of certainsensors are known and can be used as a reference.

Embodiments of the present invention that use satellite sources ofsignals provide notable benefits over existing GNSS systems, which areprone to jamming, spoofing and limitations on availability, particularlyin crowded urban areas. (Spoofing is the act of a malicious transmitterthat transmits a signal mimicking the GNSS signal, causing GNSSreceivers to output a different location and/or time from the actuallocation and/or time.) By contrast with GNSS systems, embodiments of thepresent invention are capable of exploiting signals from substantiallyany of the thousands of satellites that roam the heavens, and permitsatellite signals to be chosen for this purpose arbitrarily, on thebasis of convenience at the time of use. The set of satellite signalsthat is used may also be changed at any time if and when required, forexample if one or more of the signals are lost. The disclosedembodiments provide techniques for satellite-based timingsynchronization, localization and navigation that are more robust andless prone to jamming and spoofing than GNSS signals.

Unlike GNSS satellites, whose locations are precisely known and whosesignals are tightly controlled, the exact orbital parameters of mostsatellites in space are not known to great accuracy. The signals theyemit have different spectral characteristics and are typicallyunsynchronized with one another, and their transmission times are notknown in advance to the receivers. Embodiments of the present inventionovercome these limitations, while taking advantage of the superiorgeometrical diversity and availability of arbitrary satelliteconstellations relative to the GNSS constellations. These embodimentsprovide methods that use the radio-frequency (RF) signals received fromthe satellites at a network of terrestrial sensors to provide timingsynchronization and navigation for a sensor in the area of such anetwork and provide a superior alternative to GNSS for criticaloperations. Because the sensors can choose any desired group ofsatellites, the network is largely immune to jamming, and spoofing canbe detected easily due to the deviation of the spoofed signal from thenormal signal model. For sensor networks that also process terrestrialsignal sources (of known or unknown positions), the synchronizationaccuracy can be further improved by incorporating measurements ofsignals from terrestrial sources, such as broadcast signals or wirelessuser equipment or emitters localized by the system, in thesynchronization process.

In some embodiments, localization accuracy finer than the RF carrierwavelength can be achieved, and carrier phase measurements may then beused to refine the location and synchronization estimates still further.Furthermore, the large number of different satellite signals (as well asterrestrial sources) received by the network can be used to derive astable clock frequency in all of the sensors, thus enabling the sensorsto use longer coherent integration times. One of the benefits ofextending the coherent integration time is the ability to detect weaksignals in high-attenuation environments, such as indoors.

In another embodiment, the methods described in the present patentapplication are used to estimate the direction (azimuth and elevation)of received satellite signals or other distant sources. These estimatescan be applied in satellite orbital tracking, as well as ranging to thesatellite when Doppler measurements are also available.

System Description

FIG. 1 is a block diagram that schematically illustrates a sensornetwork 20, in accordance with an embodiment of the present invention.Sensors 22 (labeled S1, S2, S3, . . . ) are deployed at different,respective locations on a terrestrial surface 21. The sensors receive RFsignals transmitted by either terrestrial sources 24 (labeled T1, T2) orsatellite sources 26, or both. In some applications of network 20, thelocations of sensors 22, or at least of some of the sensors, are knownin advance, while the locations of sources 24 and/or 26 are not known apriori. In the description that follows, it will be assumed, for thesake of clarity and convenience, that it is the locations of sensors 22that are known and used in synchronization and source localization; butthe methods described below may be modified in a straightforward mannerto perform synchronization and localization based on known sourcepositions.

Sensors 22 are connected by a communication network 30 to asynchronization controller 28, which processes TOA data that is measuredand conveyed by the sensors over network 30 in order to performsynchronization, localization, and other related functions that aredescribed herein. Network 30 may comprise any suitable sort ofcommunication network that is known in the art. For example, network 30may comprise a dedicated, wireline or wireless network. Alternatively oradditionally, network 30 may comprise a public network, such as theInternet or a cellular communication network. Although controller 28 isshown and described here as a centralized unit, the functions of thecontroller may, additionally or alternatively, be performed in acentralized or distributed fashion by processors at different locationson network 30, including processors that are coupled to and associatedwith one or more of sensors 22. For purposes of the basic methodsdescribed below, it is sufficient that sensors 22 transmit theirrespective TOA data over network 30 to controller 28, without returntransmissions of any sort. Return transmission may be useful in someapplications, however, for conveying control, synchronization and/ornavigation information, for example, from controller 28 to sensors 22.

Although sensor network 20 and the embodiments described hereinbelowoperate on the basis of RF signals transmitted by sources 24 and 26, inalternative embodiments (not shown in the figures), sensors 22 mayreceive and measure times of arrival of signals of other types. Forexample, sources 24 may comprise acoustic emitters, and sensors 22 maycomprise acoustic receivers, which are used in localizing the emittersusing the methods described below.

FIG. 2 is a block diagram that schematically shows components of sensor22, in accordance with an embodiment of the present invention. Areceiver 32 receives and amplifies RF signals via an antenna 31 fromsources 24 and/or 26. A local processor 34 identifies the incomingsignals and measures their times of arrival, based on clock signalsprovided by a local clock 36. As noted earlier, the TOA measurements maybe simple local measurements or differential (DTOA) measurements.Processor 34 transmits the TOA measurement results for the receivedsignals from each of the relevant sources to controller 28 over network30 via a communications interface 38. Any suitable communicationprotocol, either standard or proprietary, may be used for this purpose.

As noted earlier, local clock 36 need not be particularly accurate andis generally not synchronized with the local clocks of other sensors.Any suitable sort of oscillator, as is known in the art, may thus beused in generating the clock signal according to which processor 34makes its TOA measurements. If one or more of sensors 22 has access toabsolute-time clock, then this information can be used to propagate theabsolute time among the sensors in sensor network 20 aftersynchronization is established.

FIG. 3 is a block diagram that schematically shows details ofsynchronization controller 28, in accordance with an embodiment of thepresent invention. Controller 28 receives TOA data over network 30 via acommunication interface 40. A processor 42 stores relevant data in amemory 44, and computes synchronization and localization results asdescribed hereinbelow. These results may be output via a user interface(UI) 46, for example, or transmitted to sensors 22 or to other items ofequipment via interface 40.

Processor 42 and the other components of controller 28 may be parts of ageneral-purpose computer, which is programmed in software to carry outthe functions that are described herein. This software may be downloadedto controller 28 in electronic form, over a network, for example.Alternatively or additionally, the software may be stored in a tangible,non-transitory computer-readable medium, such as optical, magnetic, orelectronic memory (which may be incorporated in memory 44). Furtheralternatively or additionally, at least some of the functions ofprocessor 42 may be carried out by dedicated or programmable hardwarelogic.

Methods of Operation

FIG. 4 is a flow chart that schematically illustrates a method for clocksynchronization and localization in sensor network 20, in accordancewith an embodiment of the present invention. This method synchronizesand localizes a system of M sensors 22 and N sources 24, 26. Certainaspects will be described below initially with reference only toterrestrial sources 24, and will then be expanded to handle satellitesources 26, as well.

To begin the procedure, the clock offset vector and cost function valuesare initialized, at an initialization step 50. In each iteration n,controller 28 estimates the location of each source t=1, 2, . . . , Nbased on the estimated clock offsets and the current value of the costfunction, at a position estimation step 52. Then, given the estimatedsource locations, controller 28 updates the clock offset estimates, at aclock update step 54.

To determine whether the estimates have converged sufficiently toterminate the computation, controller 28 compares the new cost value tothe cost value from the previous iteration, at a cost evaluation step56. Specifically, the controller may evaluate whether the change in costsince the previous iteration is less than a certain threshold. If not, nis incremented, and control returns to step 52 for a further iteration.

If the computation is found to have converged at step 56, controller 28synchronizes the clocks of sensors 22 based on the current value of theestimated clock offsets, at a synchronization step 58. Depending onapplication requirements, the controller may either send actualsynchronization values to the sensors, or it may simply record and applyappropriate corrections to the clock readings received from the sensors.The controller may also output or otherwise apply the current estimateof the source location values, at a location output step 60.

In dealing with satellite sources, it is assumed that sensors 22 aresufficiently close together (typically not more than a few tens ofkilometers apart) so that ionospheric and tropospheric effects on thesatellite signals can be considered to be uniform. Otherwise,measurements should be corrected by tropospheric and ionosphericcorrection models or measurements. If the satellite signals occupy alarge bandwidth or are scattered over a wide frequency band, theionospheric correction can be estimated accurately, as is known in theart. Alternatively or additionally, comparison of signals may be limitedto measurements from sets of sensors that are close to one another, thusallowing the entire sensor set to be deployed on a wider area.

The subsections that follow will demonstrate a number of differentimplementations of the method of FIG. 4 to synchronize the sensors:

-   -   First, methods using terrestrial sources assuming no skew at the        sensor clocks.    -   Second, methods for handling clock skew in the sensors using        terrestrial sources. (Estimating the skew using such methods is        unnecessary if satellite sources are present, as correcting the        skew using such sources is straightforward.    -   Finally methods for synchronizing the sensors using both        satellite and terrestrial sources.        Note that each subsection has its own notation.

Localization and Passive Sensor Network Synchronization UsingTerrestrial Sources—No Clock Skew

Consider M sensors and N sources. The s-th sensor has an internal clockoffset denoted by δ_(S). We use sensor number 1 as a reference sensorand therefore we set δ₁=0. Let m_(t) ^((s)) be the TOA measurementperformed by the s-th sensor for the t-th source. Denote by q_(s), p_(t)the known coordinates vector of the sensor and the unknown coordinatesvector of the source, respectively. Let τ_(t) be the unknown transmittime. Then

$\begin{matrix}{m_{t}^{(s)} = {{\frac{1}{c}{{q_{s} - p_{t}}}} + \tau_{t} + \delta_{s} + e_{t}^{(s)}}} & (1.1)\end{matrix}$where e_(t) ^(S) is a random measurement error with zero mean and astandard deviation of σ_(t) ^((s)). Here, for reasons of clarity, weassume σ_(t) ^((s))=σ_(t), or more precisely, the noise is assumed to beindependent and identically distributed (i.i.d). Under this model theproblem of self-synchronized localization can be defined as follows:Using all measurements m_(t) ^((s)) made by the sensor network estimateδ_(S), p_(t)∀s,t where s=1, 2, . . . M; t=1, 2, . . . N.

The self-synchronized localization problem, in a planar geometry, has3N+M−1 unknowns: 2N unknown source coordinates, N unknown transmit timesand M−1 sensor time offsets w.r.t. the reference sensor. If all sensorsreceive the signals of all the sources, we have MN TOA measurements.Thus, the problem is well posed if MN≥3N+M−1. Define the vectorsassociated with source tm _(t)

[m _(t) ⁽¹⁾ , m _(t) ⁽²⁾ . . . m _(t) ^((M))]^(T)r _(t)

[r _(t) ⁽¹⁾ , r _(t) ⁽²⁾ . . . r _(t) ^((M))]^(T)r _(t) ^((s))

1/c∥q _(S) −p _(t)∥1_(M)

[1, 1, . . . 1]^(T)δ

[δ₁, δ₂ . . . δ_(M)]^(T)e _(t)

[e _(t) ⁽¹⁾ , e _(t) ⁽²⁾ . . . e _(t) ^((M))]^(T).

Using the above definitions the measurements defined in (1.1) become,m _(t) =r _(t)+1_(M)τ_(t) +δ+e _(t).  (1.2)

When the error vector is Gaussian i.i.d the Maximum Likelihood costfunction is

$\begin{matrix}{Q{\sum\limits_{t = 1}^{N}{{{m_{t} - r_{t} - {1_{M}\tau_{t}} - \delta}}^{2}.}}} & (1.3)\end{matrix}$

The cost function is minimized w.r.t. the time offsets by

$\begin{matrix}{{\hat{\tau}}_{t} = {{\left( {1_{M}^{T}1_{M}} \right)^{- 1}1_{M}^{T}\left( {m_{t} - r_{t} - \delta} \right)} = {\frac{1}{M}1_{M}^{T}{\left( {m_{t} - r_{t} - \delta} \right).}}}} & (1.4)\end{matrix}$Substituting this result back in (1.3) yields

$\begin{matrix}\begin{matrix}{Q = {\sum\limits_{t = 1}^{N}{{\left( {I_{M} - {\frac{1}{M}1_{M}1_{M}^{T}}} \right)\left( {m_{t} - r_{t} - \delta} \right)}}^{2}}} \\{= {\sum\limits_{t = 1}^{N}{{B\left( {m_{t} - r_{t} - \delta} \right)}}^{2}}} \\{{= {\sum\limits_{t = 1}^{N}{\left( {{\overset{\sim}{m}}_{t} - {\overset{\sim}{r}}_{t} - \overset{\sim}{\delta}} \right)}^{2}}},}\end{matrix} & (1.5) \\{where} & \; \\{{{BI_{M}} - {\frac{1}{M}1_{M}1_{M}^{T}}}{{\overset{\sim}{m}}_{t}{Bm}_{t}}{{\overset{\sim}{r}}_{t}{Br}_{t}}{\overset{\sim}{\delta}B\;{\delta.}}} & (1.6)\end{matrix}$Collecting all vectors for t=1, 2, . . . , N we getm

[{tilde over (m)} ₁ ^(T) , m ₂ ^(T) . . . m _(N) ^(T)]^(T)r

[{tilde over (r)} ₁ ^(T) , {tilde over (r)} ₂ ^(T) . . . {tilde over(r)} _(N) ^(T)]^(T)e

[{tilde over (e)} ₁ ^(T) , {tilde over (e)} ₂ ^(T) . . . {tilde over(e)} _(N) ^(T)]^(T)C

[I _(M) , I _(M) . . . I _(M)]^(T).  (1.7)Thus, the cost function for minimization is nowQ=μm−r−C{tilde over (δ)}∥ ²  (1.8)The minimum w.r.t. {tilde over (δ)} is obtained for

$\begin{matrix}{\overset{\hat{\sim}}{\delta} = {{\left( {C^{T}C} \right)^{- 1}{C^{T}\left( {m - r} \right)}} = {\frac{1}{N}{C^{T}\left( {m - r} \right)}}}} & (1.9)\end{matrix}$Substituting back in (1.8) we get

$\begin{matrix}{{Q = {{\left\lbrack {I - {\frac{1}{N}{CC}^{T}}} \right\rbrack\left( {m - r} \right)}}^{2}},} & (1.10)\end{matrix}$

Thus, the search is “only” over {P_(t)}_(t=1) ^(N) (2N unknowns). Toensure that the number of equations is not smaller than the number ofunknowns, at least two sources and five sensors are required. Thus,solving the above equation for two sources requires a search within afour-dimensional space.

Iterative MLE

Reconsider equation (1.5) and note that if δ is given the cost functionis separable into independent components, one for each source. Thisenables us to locate each source separately, given the assumed δ vector.Therefore, the problem can be solved by successive optimization oftwo-dimensional problems. Using this fact, we can devise the followingiterative algorithm for the estimation of the clock offsets and thesources position:

1. Initialize: δ⁽⁰⁾←0; Q⁽⁻¹⁾←inf; n←0

2. Q^((n))←0

3. For each source t=1, 2, . . . , N:

-   -   (a) Estimate p_(t) given δ^((n)) using equation (1.5)    -   (b) Q^((n))←Q^((n))+∥({tilde over (m)}_(t)−{tilde over        ({circumflex over (r)})}_(t)−{tilde over (δ)}^((n)))∥²

${4.\mspace{14mu}\left. {\overset{\sim}{\delta}}^{({n + 1})}\longleftarrow\frac{1}{N} \right.{C^{T}\left( {m - r} \right)}};{{See}\mspace{14mu}{{eq}.\mspace{14mu}(1.9).}}$

5. If Q_(Δ) ^((n))=|Q^((n))−Q^((n−1))|<ϵ then stop

6. Else n←n+1, go to step 2

Note that if the position of one source is known, then the clock offsetsof the sensors that sense it can be estimated directly at the accuracyof the raw measurement. Additional refinement can then be achieved byusing that estimation as the initial guess for the δ vector. A possiblesource with known location can be a sensor that is also allowed totransmit.

The above algorithm can be used for network-wide synchronization andlocalization as long as the number of measurements is not smaller thanthe number of unknowns.

Synchronization by Convex Optimization

Define the M (M−1)/2×M difference matrix, D, where each row consists ofa single 1, a single −1 and zeros, such that Dm_(t) is the vector whoseentries are all the possible differences of the entries of m_(t). Defineδ′ such that, δ=[0 δ′^(T)]^(T) and {tilde over (D)} as the differencematrix, D, with the first column omitted. Applying D to both side ofequation (1.2) we get

$\begin{matrix}\begin{matrix}{{Dm}_{t} = {{Dr}_{t} + {D\;\delta} + {De}_{t}}} \\{{= {{Dr}_{t} + {\overset{\sim}{D}\delta^{\prime}} + {De}_{t}}},}\end{matrix} & (1.11)\end{matrix}$since D1_(M)=0, and δ₁=0 by definition. Note that Dr_(t) is the vectorholding all possible Time Difference of Arrival (TDOA) measurements.Using the triangle inequality we have |r_(t) ^((i))−r_(t)^((j))|<d^((i,j)), where d^((i,j)) is the distance between sensor i andsensor j. Collecting all possible pairs, (i, j), the triangle inequalitybecomes the component wise inequality −d

Dr_(t)

d, where d is a vector whose entries are d^((i,j)). Applying this resultto (1.11) we get,{tilde over (D)}δ′

Dm _(t) −De _(t) +d−{tilde over (D)}δ′

−Dm _(t) +De _(t) +d.  (1.12)

The set of δ′ vectors satisfying the above constraints form a polyhedronin

^(M-1) space.

Assuming that De_(t) is negligible, one can find a δ vector thatsatisfies the constraints by solving the Linear Programming (LP)feasibility problem:

Program 1.1 Clock Offset Feasibility Test

Find: δ

Subject to: ±{tilde over (D)}δ′

±Dm_(t)+d, t=1, . . . , N

This problem will return an arbitrary point in the feasible set definedby the constraints. A center point in the set would be favorable sinceits distance to the real vector will be on average less than any pointon the boundary. One way of getting a center point of the polyhedron isfinding the Maximum Volume Inscribed Ellipsoid (MVIE), which finds thecenter coordinates δ′_(C) of the biggest ellipsoid s defined as:ϵ={Gu+δ′ _(c) |∥u∥≤1, uϵ

^(M-1) , GϵS _(M-1) ⁺⁺}  (1.13)that can fit inside of the polyhedron defined by the constraints. HereS_(M-1) ⁺⁺ stands for positive definite M−1×M−1 symmetric matrix. Thenthe centering program is defined by:Program 1.2 Maximum Volume Inscribed Ellipsoid Centering on ClockOffsetsMinimize: log det G⁻¹Subject to: ∥Ga_(i) ^(T)∥+a_(i)δ′_(c)≤b_(t) _(i)

-   -   t=1, . . . , N; i=1, . . . , M(M−1)        Where:        A=[{tilde over (D)} ^(T) −{tilde over (D)} ^(T)]^(T)        b _(t)=[(Dm _(t) +d)^(T) (−Dm _(t) +d)^(T)]^(T)  (1.14)

a_(i) is the i-th row of A, and b_(t) _(i) is the i-th element of b_(t).Note that this optimization problem also returns G, which describes themaximal volume ellipsoid shape. We will use G^(T)G as a covariancematrix approximation for the estimation error of δ′.

The measure of uncertainty of the estimation of the 5 vector is usefulin a number of ways. As noted earlier, given a source of known location,the clock offsets are readily found at raw measurement accuracy.Therefore, if we have a subset of more than three clock offsetsestimates with high certainty, we can get a good source positionestimation and thus a high-certainty estimation of the clock offsets forall the sensors that received a transmission from this source.

When the noise is significant we cannot discard the term De_(t) from theconstraints in (1.12). If we collect all the constraints related to agiven sensor pair we get:δ_(i)−δ_(j)≤min({D ^((i,j)) m _(t) −D ^((i,j)) e _(t) +d ^((i,j))}_(t=1)^(N))δ_(i)−δ_(j)≥max({D ^((i,j)) m _(t) −D ^((i,j)) e _(t) −d ^((i,j))}_(t=1)^(N))where D^((i,j)) is the row in D that relates to the i, j delta pair. Themeasurement error vector, e_(t), is unknown. We can ignore themeasurement error and still use the min and max operators, but thisapproach will result in large errors.

Instead, we introduce a function that encourages small errors, thePseudo-Likelihood function, which provides the likelihood of a value ofΔ_(ij)

δ_(i)−δ_(j) given the constraints imposed on d^((i,j)) for each sensorpair using all measurements made by the pair. The Pseudo-Likelihoodfunction ƒΔ_(ij) of Δ_(ij) is defined by

M_(t)^((i, j))min²(m_(t)^((i)) − m_(t)^((j)) − (x + d^((i, j))), m_(t)^((i)) − m_(t)^((j)) − (x − d^((i, j))))$L_{t}^{({i,j})}\left\{ {\begin{matrix}1 & {{{if}\mspace{14mu}{{x - \left( {m_{t}^{(i)} - m_{t}^{(j)}} \right)}}} < d^{({i,j})}} \\{\exp\left\{ {- \frac{M_{t}^{({i,j})}}{4\sigma^{2}}} \right\}} & {otherwise}\end{matrix}{f_{\Delta_{ij}}(x)}\frac{\prod\limits_{t \in T_{ij}}L_{t}^{({i,j})}}{\int_{- \infty}^{\infty}{\prod\limits_{t \in T_{ij}}L_{t}^{({i,j})}}}} \right.$where T_(ij) is the set of sources that were measured by both sensor iand sensor j, which here include all sources.

Using the pseudo-likelihood function we can obtain thepseudo-probability that Δ_(ij)ϵ(l_(ij), h_(ij)). We select Δϵ(l_(ij),h_(ij)) so that Δ_(ij) is within the interval with pseudo-probability1−α. Since the pseudo-likelihood is normalized, it has the properties ofa probability density function, therefore:

∫_(l_(ij))^(h_(ij))f_(Δ_(ij))(x)𝕕x = 1 − αAs there is a degree of freedom in the choice of l_(ij) and h_(ij),these elements will be chosen such that ∥h_(ij)−l_(ij)∥ is minimized.

When noise is introduced, the ML estimator is such that it is notnecessary to estimate the full δ vector, but only its projection Bδ.Estimating Bδ cannot be achieved directly by using MVIE since the domainof Bδ has M−1 non-singular dimensions, making its volume in

^(M) always zero for any set of constraints. Therefore any inscribedellipsoid volume will always be zero, making any volume maximizationimpossible. To overcome this limitation, we apply the MVIE only to theprojection subspace defined by B. This subspace is the span of theeigenvectors of B corresponding to the unity eigenvalues.

We define the matrix V, whose columns are the eigenvectors of B whosecorresponding eigenvalues are different from zero. Collecting thelikelihood functions ƒΔ_(ij) for all i, j, we can derive an estimate forthe Bδ vector using MVIE centering. Estimating the Bδ vector using MVIEcentering under the constraints imposed by the pseudo-likelihood measureand its inherent singularity requires a modification of Program 1.2:

Program 1.3 Maximum Volume Inscribed Ellipsoid Centering on ClockOffsets Under Pseudo-Likelihood Imposed Bounds

Minimize: log det G⁻¹

Subject to: ∥Ga_(i) ^(T)∥+a_(i){tilde over (δ)}′_(c)≤b_(i)

-   -   i=1, . . . , M(M−1)

Return: {tilde over (δ)}_(c)←v{tilde over (δ)}′_(c)

-   -   {tilde over (G)}←VGV^(T)

Where:A=[v ^(T) D ^(T) −v ^(T) D ^(T)]^(T)b=[h _(Δ) ^(T)(α) l _(Δ) ^(T)(α)]^(T)  (1.15)

b_(i) is the i-th element of b; h_(Δ) ^(T)(α) and l_(Δ) ^(T)(α) are theupper and lower bound vectors imposed by α, via the pseudo-likelihoodfunction on all delta pairs; {tilde over (δ)}_(c) is the MVIE estimateof Bδ; and {tilde over (G)} describes the uncertainty ellipsoid of theestimation.

Note that A in Program 1.3 is not sparse. If δ′_(c) were to be estimatedusing the MVIE instead of {tilde over (δ)}′_(c), then A would be sparse,and would therefore enable lower complexity computation without loss ofaccuracy.

If Program 1.3 returns the non-feasible token, then the value of a isdecreased, and the solution is repeated.

MVIE gives an initial estimation of the target position by embedding theestimated vector in the source measurement model, which gives:Bm _(t) =B(r _(t) +e _(t))+{tilde over (δ)}_(c) +e{tilde over (δ)}_(c)  (1.16)where e_({tilde over (δ)}) _(c) is the MVIE estimation error of {tildeover (δ)}. If we model the MVIE estimation error as an independent noisewith the covariance matrix G^(T)G, then the initial source positionestimate can be found by solving the following weighted least squares(WLS) optimization problem:

$\left. {{\underset{x_{t},y_{t}}{argmin}\left( {{B\left( {m_{t} - r_{t}} \right)} - {\overset{\sim}{\delta}}_{c}} \right)}^{T}{\Sigma^{\dagger}\left( {{B\left( {m_{t} - r_{t}} \right)} - {\overset{\sim}{\delta}}_{c}} \right)}} \right)$where Σ=G^(T)G+Bσ_(t) ², and σ_(t) ² is the measurement noise varianceof the source. After locating all of the sources using the aboveestimator, the computation can jump to step 4 of the iterative MLalgorithm presented above and continue iterating until convergence isachieved. Two other possible ways to proceed are

-   -   1. Perform gradient ascent on the full ML cost using the MVIE as        the initial guess;    -   2. Use the MVIE initial estimate as the a priori information in        maximum a posteriori estimation (MAP).

Using G^(T)G as a covariance matrix is conservative, since it describesthe ellipsoid within which the clock offsets vector almost always lie.Therefore adding a constant factor to the matrix to convert it from acertainty ellipsoid to a covariance matrix is appropriate. Defining theconstant factor as κ>1, the covariance matrix for the MVIE clock offsetestimation is

$\frac{G^{T}G}{\kappa}.$Methods are known in the art for computing κ to convert certaintyellipsoid to a covariance matrix of a multi-variate normal distribution,but κ can also be determined empirically.

The MVIE/ML estimator performs well under the imposed constraints. TheMVIE/ML estimator even after only one ML iteration comes very close tothe converged ML performance. This feature makes the combined MVIE/MLestimator very attractive in terms of processing time. Comparing theresults of the combined MVIE/ML estimator to the known clock offsetcase, we observe that even in the first iteration, the penalty of notknowing the clock offsets is not very large. Moreover, even when σ_(n)is on the order of the distance between the sensors (which is a highlynonlinear configuration), the MVIE/ML estimator converges successfully.Moreover the MVIE can also be used as a priori information in a MAPestimator.

To simplify the description above, it was assumed that all sensorsintercept all the sources. It is straightforward to extend the aboveresults to the more general case.

Note also that the ML estimator cannot use measurements of sources thatare received by three sensors or less, but MVIE can.

One further point to note is that clock bias is actually the measurementbias, meaning that it contains any propagation delays inside of thesensor. For example, a sensor that has an antenna connected with a longcable will include the cable propagation time as part of the bias. Insensing systems, this effect should generally be calibrated. The methoddescribed above provides sensor self-calibration of such delays.

2D Bounds on DTOA Measurements

In the previous section, the MVIE used constraints based on measurementsof a source by a pair of sensors. This section will present constraintsbased on measurements of a source by a sensor triplet. Methods toincorporate these constraints into the MVIE will also be introduced.

Let there be three sensors with position coordinates q_(i). Setdirection vectors v_(ij) as:v _(ij)

q _(j) −q _(i)  (2.1)Define two angles:α

∠(v ₀₁ ,v ₀₂)β

∠(v ₁₂ ,−v ₀₁)  (2.2)where the convention is that for vectors μ, ν the angle ∠(μ, ν) ispositive if the ν direction is restored by counterclockwise rotation ofμ by less than a π rad turn. The noiseless DTOA measurement made by apair sensors on a target positioned at p is defined, in metric units,as:ψ_(ij)(p)=∥p−q _(j) ∥−∥p−q _(i)∥  (2.3)

Measurement space is defined as the R² space (ψ₀₁, ψ₀₂) A nine-segment{u_(i)}_(i=0) ⁸ path in space is established to map to the bound inmeasurement space. The path will traverse the plane spanned by vectorsv₀₁, v₀₂, as follows:u ₀(t ₀)=q ₀ −t ₀ v ₀₁u ₁(t ₁)=q ₀ −t ₁ v ₀₂u ₂(t ₂)=q ₁ +t ₂ v ₀₁u ₃(t ₃)=q ₁ −t ₃ v ₁₂u ₄(t ₄)=q ₂ +t ₄ v ₀₂u ₅(t ₅)=q ₂ +t ₅ v ₁₂u ₆(φ₀)=q ₀ +rT(φ₀)v ₀₁, φ₀ϵ[0, α]u ₇(φ₁)=q ₀ +rT(φ₁)v ₀₁, φ₁ϵ[π−β, π]u ₈(φ₂)=q ₀ +rT(φ₂)v ₀₁, φ₂ϵ[π+α, 2π−β]t _(i)>0, ∀i  (2.4)where r goes to infinity, and T(φ) is the rotation matrix for a counterclock wise rotation of φ degrees.

Traveling the path segments in space translates to the following path inmeasurement space:h ₀(l ₀)=l ₀ [|v ₀₁ |, |v ₀₂|]^(T)+(1−l ₀)[|v ₀₁ |, |v ₀₂|cos(α−π)]^(T)h ₁(l ₁)=l ₁ [|v ₀₁ |, |v ₀₂|]^(T)+(1−l ₁)[|v ₀₁|cos(π+α), |v ₀₂|]^(T)h ₂(l ₂)=l ₂ [−|v ₀₁ |, |v ₁₂ |−|v ₀₁|]^(T)+(1−l ₂)[−|v ₀₁ |, −|v ₀₂|cosα]h ₃(l ₃)=l ₃ [−|v ₀₁ |, |v ₁₂ |−|v ₀₁|]^(T)+(1−l ₃)[−|v ₀₁|cos β, −|v₀₂|cos(α+β)]^(T)h ₄(l ₄)=l ₄ [|v ₁₂ |−|v ₀₂ |, −|v ₀₂|]^(T)+(1−l ₄)[−|v ₀₁|cos α, −|v₀₂|]^(T)h ₅(l ₅)=l ₅ [|v ₁₂ |−|v ₀₂ |, −|v ₀₂|]^(T)+(1−l ₅)[|v ₀₁|cos β, |v₀₂|cos(α+β)]^(T)l _(i)ϵ[0,1]∀ih ₆(φ₀)=−[|v ₀₁|cos φ₀ , |v ₀₂|cos(α−φ₀)]^(T)φ₀ϵ[0,α]h ₇(φ₁)=−[|v ₀₁|cos φ₁ , |v ₀₂|cos(α−φ₁)]^(T)φ₁ϵ[π−β, π]h ₈(φ₂)=−[|v ₀₁|cos φ₂ , |v ₀₂|cos(α−φ₂)]^(T)φ₂ ϵ[π+α,2π−β]  (2.5)

The elementary bound, derived from the triangle inequality,−|ν_(ij)|≤ψ_(ij)≤|ν_(ij)|, is directly responsible for the linear partof the bounding path. Therefore the additional information gained byconsidering two-dimensional bounds is found in the nonlinear part of thebounding path (namely h_(6,7,8)).

To avoid the complexity of incorporating the nonlinear bound in theconvex optimization described above, it is possible to find the tangentsto the nonlinear part of the bound. Each tangent will act as ahalf-space bound on the offset values. The more tangents we use, thetighter the bound will be, but at the cost of greater complexity.

Bound Tangents

The nonlinear part of the bound is composed of three sections of theellipse:ψ₀₂=−|ν₀₂|cos(α−φ)ψ₀₁=−ν₀₁|cos(φ)φϵ[0, α]∪[π−β, π]∪[+α, 2π−β]  (2.6)Taking the partial derivative of ψ₀₁ and ψ₀₂ by φ gives:

$\begin{matrix}{{\frac{\partial\psi_{02}}{\partial\varphi} = {{- {v_{02}}}{\sin\left( {\alpha - \varphi} \right)}}}{\frac{\partial\psi_{01}}{\partial\varphi} = {{v_{01}}{\sin(\varphi)}}}} & (2.7)\end{matrix}$Dividing the two expressions results in:

$\begin{matrix}{{\frac{\partial\psi_{02}}{\partial\psi_{01}}❘_{\varphi}} = {- \frac{{v_{02}}{\sin\left( {\alpha - \varphi} \right)}}{{v_{01}}{\sin(\varphi)}}}} & (2.8)\end{matrix}$The tangent equation is then:

$\begin{matrix}\begin{matrix}{\psi_{02} = {\psi_{02}❘_{\varphi}{{{+ \left( {{\psi_{01} - \psi_{01}}❘_{\varphi}} \right)}\frac{\partial\psi_{02}}{\partial\psi_{01}}}❘_{\varphi}}}} \\{= {{{- {v_{02}}}\left( {{\cos\left( {\alpha - \varphi} \right)} + \frac{\sin\left( {\alpha - \varphi} \right)}{\tan\;\varphi}} \right)} - {\frac{{v_{02}}{\sin\left( {\alpha - \varphi} \right)}}{{v_{01}}{\sin(\varphi)}}\psi_{01}}}} \\{{q_{0}(\varphi)} + {{q_{1}(\varphi)}{\psi_{01}.}}}\end{matrix} & (2.9)\end{matrix}$

As noted above, there are three segments of the ellipse that are part ofthe bound. We take K tangents for each segment, evenly spaced in termsof φ so that:

$\begin{matrix}{{h_{6}:{\varphi \in \left\{ \frac{\alpha\; k}{K + 1} \right\}_{k = 1}^{K}}}{h_{7}:{\varphi \in \left\{ {\pi - \beta + \frac{\beta\; k}{K + 1}} \right\}_{k = 1}^{K}}}{h_{8}:{\varphi \in \left\{ {\pi + \alpha + \frac{\left( {\pi - \beta - \alpha} \right)k}{K + 1}} \right\}_{k = 1}^{K}}}} & (2.10)\end{matrix}$

The segments h_(6,7) are always in the lower part of the bound, so thattheir tangents pose a lower bound, while h₈ is always on top, so thatits tangents pose an upper bound.

One-Sided Pseudo Likelihood and Tangent Constraints on the Clock Offsets

The bounds on a pair of DTOA measurements made on a single source arerelated to the set of possible clock offsets. Let the measurement madeon a source t by sensor s be:m _(t) ^((s)=r) _(t) ^((s))+τ_(t)+δ_(s) +e _(t) ^((s))  (2.11)

Define {tilde over (m)}^((k,s,p))

m_(t) ^((p))−m_(t) ^((k))−q₁(φ)(m_(t) ^((s))−m_(t) ^((k))) andΔ_(ij)=δ_(j)−δ_(i). Expanding this expression gives:

$\begin{matrix}\begin{matrix}{{\overset{\sim}{m}}_{t}^{({k,s,p,\varphi})}\overset{\Delta}{=}{m_{t}^{(p)} - m_{t}^{(k)} - {{q_{1}(\varphi)}\left( {m_{t}^{(s)} - m_{t}^{(k)}} \right)}}} \\{{= {\psi_{t}^{({pk})} + \Delta_{kp} - {{q_{1}(\varphi)}\left( {\psi_{t}^{({sk})} + \Delta_{sk}} \right)} + {\overset{\sim}{e}}_{t}}}\;} \\{= {\psi_{t}^{({pk})} - {{q_{1}(\varphi)}\psi_{t}^{({sk})}} + \Delta_{kp} - {{q_{1}(\varphi)}\Delta_{sk}} + {\overset{\sim}{e}}_{t}}} \\{\underset{<}{\overset{>}{=}}{{q_{0}(\varphi)} + \Delta_{kp} - {{q_{1}(\varphi)}\Delta_{sk}} + {\overset{\sim}{e}}_{t}}} \\{\overset{\Delta}{=}{{q_{0}(\varphi)} + {\overset{\sim}{\Delta}}_{ksp}^{(\varphi)} + {\overset{\sim}{e}}_{t}}}\end{matrix} & (2.12)\end{matrix}$where {tilde over (e)}_(t)

e_(t) ^((p))+(q₁(φ)−1)e_(t) ^((k))−q₁(φ)e_(t) ^((s)) and ψ_(t)^((ks))=r_(t) ^((s))−r_(t) ^((k)). The last equation defines therelation between the bound on the DTOA measurements and the clockoffsets. The direction of the inequality is determined by the value ofφ: For φϵ[0, α]∪[π−β, π] it will be ≥ and for φϵ[π+α, 2π−β] it will be≤.

To fuse all of the measurements made by all the sensors we again turn tothe pseudo likelihood function. The pseudo likelihood functionƒ_({tilde over (Δ)}) _(ksp) ^((φ)) of {tilde over (Δ)}_(ksp) ^((φ)) isdefined by:

${\overset{\sim}{m}}_{t}^{({k,s,p,\varphi})}\overset{\Delta}{=}{m_{t}^{(p)} - m_{t}^{(k)} - {{q_{1}(\varphi)}\left( {m_{t}^{(s)} - m_{t}^{(k)}} \right)}}$${\Psi(x)}\overset{\Delta}{=}{\exp\left\{ {- \frac{\left( {{\overset{\sim}{m}}_{t}^{({k,s,p,\varphi})} - {q_{0}(\varphi)}} \right)^{2}}{2{\overset{\sim}{\sigma}}_{t}^{2}}} \right\}}$$L_{{\overset{\sim}{m}}_{t}^{({k,s,p,\varphi})}}\overset{\Delta}{=}\left\{ {{\begin{matrix}1 & {{{if}\mspace{14mu} x}\underset{<}{\overset{>}{=}}{{\overset{\sim}{m}}_{t}^{({k,s,p,\varphi})} - {q_{0}(\varphi)}}} \\{\Psi(x)} & {otherwise}\end{matrix}{f_{{\overset{\sim}{\Delta}}_{k,s,p,\varphi}}(x)}}\overset{\Delta}{=}\frac{\prod\limits_{t \in T_{ksp}}\; L_{{\overset{\sim}{m}}_{t}^{({k,s,p,\varphi})}}}{\int_{- \infty}^{\infty}{\prod\limits_{t \in T_{ksp}}\; L_{{\overset{\sim}{m}}_{t}^{({k,s,p,\varphi})}}}}} \right.$where {tilde over (σ)}_(t) ²=2σ_(t) ² (1+q₁(φ)), and T_(ksp) is the setof all sensor triplets. The direction of the inequality is determined bythe value of φ as before. The pseudo likelihood function yields theconstraint

${\overset{\sim}{\Delta}}_{ksp}^{(\varphi)}\overset{>}{\underset{<}{=}}h_{ksp}^{\varphi}$w.p. 1−ϵ. This is a linear constraint on the clock offsets that can beincorporated into the MVIE optimization problem.

Since this is a one-sided constraint, some a priori knowledge is neededon the bound of {tilde over (Δ)}_(ksp) ^((φ)) on its other side. Forthis purpose, the bounds can be derived from the elementary bounds, orthe lowest feasible value can be used.

Localization and Passive Sensor Network Synchronization Under a LinearClock Model

In this section, the above methods are expanded to include clock skewestimation by convex methods. With slight approximations, the skewestimation is independent of the offset estimation and can be achievedby similar methods.

Consider M sensors and N sources. The s-th sensor has an internal clockwith offset denoted by δ_(s) and skew denoted by θ_(s). We use sensornumber 0 as a reference sensor, and therefore we set δ₀=0 and θ₀=1. Letm_(t) ^((s)) be the TOA measurement performed by the s-th sensor for thet-th source. Denote by q_(s), p_(t) the known coordinate vector of thesensor and the unknown coordinate vector of the source, respectively.Then r_(t) ^((s))

1/c∥q_(s)−p_(t)∥ is the propagation time between source t and sensor s.Let τ_(t) be the unknown transmit time. Thenm _(t) ^((s))=θ_(s)(r _(t) ^((s))+τ_(t))+δ_(s) +e _(t) ^((s))  (3.1)where e_(t) ^((s)) is a random measurement error with zero mean and astandard deviation of σ_(t) ^((s)).

It can be reasonably assumed that all clock skews can be described by aslight deviation from unity that is denoted by η:θ_(s)=1+η_(s)  (3.2)With reasonable values of η, it can be shown thatθ_(S) r _(t) ^((s)) ≈r _(t) ^((s))  (3.3)is a good approximation since η_(s)r_(t) ^((s))<<σ_(t) ^((s)). For thesame reason we will assume that multiplying the noise by any clock skewdoes not change the noise statistics.

Given this model, the problem of self-synchronized localization can bedefined as follows: Using all measurements m_(t) ^((s)) made by thesensor network, estimate δ_(s), θ_(s), p_(t)∀s,t where s=1, 2, . . .M−1; t=1, 2, . . . N.

The self-synchronized localization problem, in a planar geometry, has3N+2(M−1) unknowns: 2N unknown source coordinates, N unknown transmittimes, M−1 sensor time offsets, and M−1 clock skews w.r.t. the referencesensor. If all sensors receive the signals of all the sources, there areMN TOA measurements. Thus, the problem is well posed if MN≥3N+2(M−1).Define the vectors associated with source t,m _(t)

[m _(t) ⁽¹⁾ , m _(t) ⁽²⁾ . . . m _(t) ^((M-1))]^(T) ; m _(t0)

[m _(t) ⁽⁰⁾, 0 . . . 0]^(T);r _(t)

[r _(t) ⁽¹⁾ , r _(t) ⁽²⁾ . . . r _(t) ^((M-1))]^(T) ; r _(t) ^((s))

1/c∥q _(s) −p _(t)∥;1_(M)

[1, 1, . . . 1]^(T); δ

[δ₁, δ₂ . . . δ_(M-1)]^(T);θ

[θ₁, θ₂ . . . θ_(M-1)]^(T) ; e _(t)

[e _(t) ⁽⁰⁾ , e _(t) ⁽¹⁾ . . . e _(t) ^((M-1))]^(T){tilde over (δ)}

δ⊕θ; {tilde over (θ)}

1_(M-1)⊕θA _({tilde over (δ)})

[0, I _(M-1)]^(T) ; A _({tilde over (θ)}) _(t)

[0, diag[m _(t)]]^(T)where ⊕ is defined as the element-wise division operator. Using theabove definitions, the measurements defined in (3.1) become,m _(t0) =r _(t)+1Mτ _(t) +A _({tilde over (δ)}) {tilde over (δ)}−A_({tilde over (θ)}) _(t) {tilde over (θ)}+e _(t).  (3.4)

When the error vector is Gaussian, the Maximum Likelihood cost functionis

$\begin{matrix}{Q\overset{\Delta}{=}{\sum\limits_{t = 1}^{N}\;{{{m_{t\; 0} - r_{t} - {1_{M}\tau_{t}} - {A_{\overset{\sim}{\delta}}\overset{\sim}{\delta}} + {{A_{\overset{\sim}{\theta}}}_{t}\overset{\sim}{\theta}}}}^{2}.}}} & (3.5)\end{matrix}$The cost function is minimized w.r.t. the time offsets by:

$\begin{matrix}\begin{matrix}{{\hat{\tau}}_{t} = {\left( {1_{M}^{T}1_{M}} \right)^{- 1}1_{M}^{T}\left( {m_{t\; 0} - r_{t} - {A_{\overset{\sim}{\delta}}\overset{\sim}{\delta}} + {{A_{\overset{\sim}{\theta}}}_{t}\overset{\sim}{\theta}}} \right)}} \\{= {\frac{1}{M}1_{M}^{T}\left( {m_{t\; 0} - r_{t} - {A_{\overset{\sim}{\delta}}\overset{\sim}{\delta}} + {{A_{\overset{\sim}{\theta}}}_{t}\overset{\sim}{\theta}}} \right)}}\end{matrix} & (3.6)\end{matrix}$Substituting this result back in (3.5) yields:

$\begin{matrix}{\begin{matrix}{Q = {\sum\limits_{t = 1}^{N}\;{{\left( {I_{M} - {\frac{1}{M}1_{M}1_{M}^{T}}} \right)\left( {m_{t\; 0} - r_{t} - {A_{\overset{\sim}{\delta}}\overset{\sim}{\delta}} + {{A_{\overset{\sim}{\theta}}}_{t}\overset{\sim}{\theta}}} \right)}}^{2}}} \\{= {\sum\limits_{t = 1}^{N}\;{{B\left( {m_{t\; 0} - r_{t} - {A_{\overset{\sim}{\delta}}\overset{\sim}{\delta}} + {{A_{\overset{\sim}{\theta}}}_{t}\overset{\sim}{\theta}}} \right)}}^{2}}} \\{{= {\sum\limits_{t = 1}^{N}\;{\left( {{\overset{\sim}{m}}_{t} - {\overset{\sim}{r}}_{t} - {A_{t}\overset{\sim}{\gamma}}} \right)}^{2}}},}\end{matrix}{where}} & (3.7) \\{{{B\overset{\Delta}{=}{I_{M} - {\frac{1}{M}1_{M}1_{M}^{T}}}};{A_{t}\overset{\Delta}{=}\left\lbrack {{BA}_{\overset{\sim}{\delta}},{- {BA}_{{\overset{\sim}{\theta}}_{t}}}} \right\rbrack}}{{{\overset{\sim}{m}}_{t}\overset{\Delta}{=}{Bm}_{t\; 0}};{{\overset{\sim}{r}}_{t}\overset{\Delta}{=}{Br}_{t}};{\overset{\sim}{\gamma}\overset{\Delta}{=}{\left\lbrack {\overset{\sim}{\delta},\overset{\sim}{\theta}} \right\rbrack.}}}} & (3.8)\end{matrix}$Collecting all vectors for t=1, 2, . . . ,N gives:m

[{tilde over (m)} ₁ ^(T) , {tilde over (m)} ₂ ^(T) . . . {tilde over(m)} _(N) ^(T)]^(T) ; r

[{tilde over (r)} ₁ ^(T) , {tilde over (r)} ₂ ^(T) . . . {tilde over(r)} _(N) ^(T)]^(T);e

[{tilde over (e)} ₁ ^(T) , {tilde over (e)} ₂ ^(T) . . . {tilde over(e)} _(N) ^(T)]^(T) ; C

[A ₁ , A ₂ . . . A _(N)]^(T).  (3.9)

Thus, the cost function for minimization is now:Q=∥m−r−C{tilde over (γ)}∥ ²  (3.10)The minimum w.r.t. {tilde over (γ)} is obtained for{tilde over ({circumflex over (γ)})}=(C ^(T) C)⁻¹ C ^(T)(m−r)  (3.11)Substituting back in (3.10) gives:

$\begin{matrix}{{Q = {{\left\lbrack {I - {\frac{1}{N}{CC}^{T}}} \right\rbrack\left( {m - r} \right)}}^{2}},} & (3.12)\end{matrix}$

As a result, the search is over only {p_(t)}_(t=1) ^(N) (2N unknowns).To ensure that the number of equations is not smaller than the number ofunknowns, at least three sources and seven sensors are required (for theminimal number of sources). Thus, solving the above equation for threesources requires a search within a 6 dimensional space.

Iterative MLE

Reconsider equation (3.7) and note that if {tilde over (γ)} is given,the cost function is separable into independent components, one for eachsource. This separation enables controller 28 to locate each sourceseparately, given the assumed {tilde over (γ)} vector. Therefore, theproblem can be solved by successive optimization of two-dimensionalproblems. Using this fact, the following iterative algorithm can beapplied for the estimation of the clock offsets and the sourcepositions:

1. Initialize: {tilde over (γ)}⁽⁰⁾←[0_(M-1) ^(T), 1_(M-1) ^(T)]^(T);Q⁽⁻¹⁾←inf; n←0

2. Q^((n))←0

3. For each source t=1, 2, . . . , N:

-   -   (a) Estimate p_(t) given γ^((n)) using equation (3.7).    -   (b) Q^((n))←Q^((n))+∥({tilde over (m)}_(t)−{tilde over        ({circumflex over (r)})}_(t)−{tilde over (γ)}^((n)))∥²

4. {tilde over (γ)}^((n+1))←(C^(T)C)⁻¹C^(T) (m−r); See eq. (3.11).

5. If Q_(Δ) ^((n))=|Q^((n))−Q^((n−1))|<ϵ then stop

6. Else n←n+1, go to step 2

Self-Synchronized Localization by Convex Optimization

Controller 28 starts by collecting all measurements made by a pair ofsensors k, s. For a source t the measurement of sensor k is given by Eq.3.1. Taking the difference between the measurements made by sensors kand s on source t gives:m _(t) ^((k)) −m _(t) ^((s))=θ_(k)(r _(t) ^((k))+τ_(t))+δ_(k) +e _(t)^((k))−θ_(s)(r _(t) ^((s))+τ_(t))−δ_(s) −e _(t) ^((s))≈r_(t) ^((k)) −r_(t) ^((s))+τ_(t)(θ_(k)−θ_(s))+δ_(k) +e _(t) ^((k)−δ) _(s) −e _(t)^((s))=φ_(t) _(ks) +τ_(t)(θ_(k)−θ_(s))+δ_(k)−δ_(s) +e _(t) ^((k)) −e_(t) ^((s))  (3.25)where φ_(t) _(ks)

r_(t) ^((k))−r_(t) ^((s)) is the DTOA measurement on source t by sensorsk,s.

We define:Δ_(ks)

δ_(k)−δ_(s)dm _(t) _(ks)

m _(t) ^((k)) −m _(t) ^((s))  (3.26)The unsynchronized DTOA measurement for source t by sensors k and s isthen:dm _(t) _(ks) ≈φ_(t) _(ks) +τ_(t)(θ_(k)−θ_(s))+Δ_(ks) +e _(t) ^((k)) −e_(t) ^((s))  (3.27)

Taking the difference between two unsynchronized DTOA measurements ofdifferent sources gives:dm _(t) _(ks) −dm _(p) _(ks) ≈φ_(t) _(ks) −φ_(p) _(ks)+(τ_(t)−τ_(p))(θ_(k)−θ_(s))+e _(t) ^((k)) −e _(t) ^((s)) +e _(p) ^((k))−e _(p) ^((s))=φ_(t) _(ks) −φ_(p) _(ks)+(τ_(t)−τ_(p))(θ_(k)−θ_(s))+{tilde over (e)} _(t,9) ^((k,s))where {tilde over (e)}_(t,p) ^((k,s))

^(e) _(t) ^((k))−e_(t) ^((s))+e_(p) ^((k))−e_(p) ^((s)). Whenconsidering a sensor pair, the relative clock skew may be defined asΘ_(ks)

θ_(k)/θ_(s). Applying this concept gives:dm _(t) _(ks) −dm _(p) _(ks) ≈φ_(t) _(ks) −φ_(p) _(ks)+(τ_(t)−τ_(p))(Θ_(ks)−1)+{tilde over (e)} _(t,p) ^((k,s))=φ_(t) _(ks)−φ_(p) _(ks) +(τ_(t)−τ_(p))(η_(ks))+{tilde over (e)} _(t,p)^((k,s))  (3.28)

Since in Eq. 3.28, τ_(t)−τ_(p) is multiplied by the clock skew deviationη, only a crude estimation of τ_(t)−τ_(p) is needed to get an accurateclock skew estimation. The following approximation can be used for theTx time difference:

$\begin{matrix}{= {{\frac{1}{M}{\sum\limits_{k = 0}^{M - 1}\;\left\lbrack {m_{t}^{(k)} - m_{p}^{(k)}} \right\rbrack}}\overset{\Delta}{=}{\Delta\; T_{t,p}}}} & (3.29)\end{matrix}$The summation is over the common sensing sensors of t and p.

Dividing by the estimated time difference gives:

$\begin{matrix}{\frac{{dm}_{t_{ks}} - {dm}_{p_{ks}}}{\Delta\; T_{t,p}} \approx {\frac{\phi_{t_{ks}} - \phi_{p_{ks}}}{\Delta\; T_{t,p}} + \eta_{ks} + \frac{{\overset{\sim}{e}}_{t,p}^{({k,s})}}{\Delta\; T_{t,p}}}} & (3.30)\end{matrix}$Finally we get:

$\begin{matrix}{{\frac{{dm}_{t_{ks}} - {dm}_{p_{ks}}}{\Delta\; T_{t,p}} - \frac{\phi_{t_{ks}} - \phi_{p_{ks}}}{\Delta\; T_{t,p}}} \approx {\eta_{ks} + \frac{{\overset{\sim}{e}}_{t,p}^{({k,s})}}{\Delta\; T_{t,p}}}} & (3.31)\end{matrix}$The triangle inequality gives bounds on η_(ks):

$\begin{matrix}{{{\eta_{ks} + {\overset{\sim}{n}}_{t,p,k,s}} \geq {\frac{{dm}_{t_{ks}} - {dm}_{p_{ks}}}{\Delta\; T_{t,p}} - \frac{2\; d_{ks}}{c\;\Delta\; T_{t,p}}}}{{\eta_{ks} + {\overset{\sim}{n}}_{t,p,k,s}} \leq {\frac{{dm}_{t_{ks}} - {dm}_{p_{ks}}}{\Delta\; T_{t,p}} + \frac{2\; d_{ks}}{c\;\Delta\; T_{t,p}}}}} & (3.32)\end{matrix}$where d_(ks)

∥q_(s)−q_(k).

The pseudo likelihood function provides the likelihood of a value ofη_(ks) given the constraints imposed by Eq. 3.32 for each sensor pairusing all measurements made by the pair. Assuming all measurement noiseis i.i.d and has variance of σ_(t), the pseudo likelihood function ƒ_(η)_(ks) of η_(ks) is defined by:

$\begin{matrix}{{{ddm}_{t,p}^{({k,s})}\overset{\Delta}{=}{{dm}_{t_{ks}} - {dm}_{p_{ks}}}}{M_{t,p}^{({k,s})}\overset{\Delta}{=}{\min^{2}\begin{pmatrix}{{{\frac{{d{dm}}_{t,p}^{({k,s})}}{\Delta\; T_{t,p}} + \frac{2\; d_{k,s}}{c\;\Delta\; T_{t,p}} - \eta_{ks}}},} \\{{\frac{{d{dm}}_{t,p}^{({k,s})}}{\Delta\; T_{t,p}} - \frac{2\; d_{k,s}}{c\;\Delta\; T_{t,p}} - \eta_{ks}}}\end{pmatrix}}}{{\Psi\left( \eta_{ks} \right)}\overset{\Delta}{=}{\exp\left\{ {- \frac{\Delta\; T_{t,p}^{2}M_{t,p}^{({k,s})}}{8\sigma_{t}^{2}}} \right\}}}{L_{{ddm}_{t,p}}^{({k,s})}\overset{\Delta}{=}\left\{ {{\begin{matrix}1 & {{{if}\mspace{14mu}{{\eta_{ks} - \frac{{d{dm}}_{t,p}^{({k,s})}}{\Delta\; T_{t,p}}}}} < \frac{2\; d_{k,s}}{c\;\Delta\; T_{t,p}}} \\{\Psi\left( \eta_{ks} \right)} & {otherwise}\end{matrix}{f_{\eta_{ks}}(x)}}\overset{\Delta}{=}\frac{\prod\limits_{{\{{t,p}\}} \in T_{k,s}}\; L_{{ddm}_{t,p}}^{({k,s})}}{\int_{- \infty}^{\infty}{\prod\limits_{{\{{t,p}\}} \in T_{k,s}}\; L_{{ddm}_{t,p}}^{({k,s})}}}} \right.}} & (3.33)\end{matrix}$Here T_(k,s) is the set of source pairs that were measured by bothsensor k and sensor s, which is assumed here to include all sources. Itis assumed that ddm_(t,p) ^((k,s)) are uncorrelated with one another,but this approximation can easily be eliminated. From the pseudolikelihood function, we obtain E [η_(ks)], var[η_(ks)], andη_(ks)ϵ(l_(ks), h_(ks)) w.p. 1−α.

The estimation for the relative clock skew is given by {circumflex over(Θ)}_(ks)=1+{circumflex over (η)}_(ks). Therefore, bounds that werefound on η_(ks) can be translated into bounds on Θ_(ks). Ifη_(ks)ϵ(l_(ks), h_(ks)) w.p. 1−α, then Θ_(ks)ϵ(1+l_(ks), 1+h_(ks)) w.p.1−α. Taking the natural logarithm of the inequalities imposed on Θ_(ks)gives the following set of linear constraints on ln θ_(k):ln Θ_(ks)=ln θ_(k)−ln θ_(s)ϵ[ln(1+l _(ks)), ln(1+h _(ks))] w.p.1−α.  (3.34)

Collecting the linear constraints for all k, s sensor pairs, MaximumVolume Inscribed Ellipsoid (MVIE) centering can be used to find theestimate

, and the related error covariance matrix can be approximated byG_(ln θ)G_(ln θ) ^(T).

Program 3.1 Maximum Volume Inscribed Ellipsoid Centering on Clock SkewsUnder Pseudo Likelihood Imposed Bounds

Minimize: log det G_(ln θ) ⁻¹

Subject to: ∥G_(ln θ)a_(i) ^(T)∥+a_(i) ln θ_(c)≤b_(i)

-   -   i=1, . . . , M(M−1)

Where:A=[{tilde over (D)} ^(T) −{tilde over (D)} ^(T)]^(T)b=[h _(ln Θ) ^(T)(α) l _(ln Θ) ^(T)(α)]^(T)  (3.35)b_(i) is the i-th element of b, and h_(ln Θ) ^(T)(α) and l_(ln Θ)^(T)(α) are the upper- and lower-bound vectors imposed by α, via thepseudo likelihood function, on all relative clock skews.

If the program returns the non-feasible token, the value of α can bedecreased, and the solution repeated.

Note that G_(θ)=G_(ln θ), since the error is much smaller then unity,and Taylor approximation holds firmly. Therefore, the estimator of the θvector is:{circumflex over (θ)}=exp(

); G _(θ) =G _(ln θ)  (3.36)

After establishing an estimation of the clock skew vector {circumflexover (θ)} and deriving an approximation for its error covariance matrix:G_(θ)G_(θ) ^(T), the next step is to estimate the clock offset vector δ.For this purpose, the measurements made by the sensors are normalized bytheir respective clock skew estimates as follows:

$\begin{matrix}{{\overset{\sim}{m}}_{t}^{(s)} = {\frac{m_{t}^{(s)}}{{\hat{\theta}}_{s}} = {{\frac{\theta_{s}}{{\hat{\theta}}_{s}}\left( {r_{t}^{(s)} + \tau_{t}} \right)} + \frac{\delta_{s}}{{\hat{\theta}}_{s}} + e_{t}^{(s)}}}} & (3.37)\end{matrix}$

Applying Taylor approximation to

$\frac{\theta_{s}}{{\hat{\theta}}_{s}}$gives

${\frac{\theta_{s}}{{\hat{\theta}}_{s}} \approx {1 - e_{\theta_{s}}}},$where e_(θ) _(s) is the estimation error of {circumflex over (θ)}_(s).Embedding this approximation into Eq. 3.37 results in:

$\begin{matrix}\begin{matrix}{{\overset{\sim}{m}}_{t}^{(s)} = {{\left( {1 - e_{\theta_{s}}} \right)\left( {r_{t}^{(s)} + \tau_{t}} \right)} + \frac{\delta_{s}}{{\hat{\theta}}_{s}} + e_{t}^{(s)}}} \\{= {r_{t}^{(s)} + \tau_{t} + \frac{\delta_{s}}{{\hat{\theta}}_{s}} - {e_{\theta_{s}}\left( {r_{t}^{(s)} + \tau_{t}} \right)} + e_{t}^{(s)}}} \\{\approx {r_{t}^{(s)} + \tau_{t} + \frac{\delta_{s}}{{\hat{\theta}}_{s}} - {e_{\theta_{s}}\tau_{t}} + e_{t}^{(s)}}}\end{matrix} & (3.38)\end{matrix}$

The above equation states that the clock skew uncertainty is amplifiedby the source transmit (Tx) time. The Tx time can be referenced to aselected point in time. To minimize the noise amplification effect, thereference time can be set to the mean of the Tx times of all of thesources. Then τ_(t) may be estimated by:

$\begin{matrix}{{\tau_{0}\overset{\Delta}{=}{\frac{1}{N}{\sum\limits_{k = 1}^{N}\; m_{k}^{(0)}}}}{{\hat{\tau}}_{t} = {m_{t}^{(0)} - \tau_{0}}}} & (3.39)\end{matrix}$

Using the above first-order approximation of τ_(t), the equivalentmeasurement difference noise for sensors k and s is:

$\begin{matrix}{\sigma_{ks} = \sqrt{\frac{{\hat{\tau}}_{t}^{2}v_{ks}^{T}G_{\theta}G_{\theta}v_{ks}}{4} + {2\sigma_{t}^{2}}}} & (3.40)\end{matrix}$where v_(ks) ^(T) is a difference vector holding 1 and −1 in the k−1,s−1 rows respectively, and containing zeros elsewhere. (If s=1, then the−1 is missing from the vector.)Given σ_(ks) the pseudo likelihood measure can be applied to give boundson

${\overset{\sim}{\Delta}}_{ks}\overset{\Delta}{=}{{\frac{\delta_{k}}{{\hat{\theta}}_{k}} - \frac{\delta_{s}}{{\hat{\theta}}_{s}}}\overset{\Delta}{=}{{\overset{\sim}{\delta}}_{k} - {\overset{\sim}{\delta}}_{s}}}$for every s,k. The {tilde over (Δ)} pseudo likelihood measure is definedas follows:

${d{\overset{\sim}{m}}_{t}^{({i,j})}}\overset{\Delta}{=}{{\overset{\sim}{m}}_{t}^{(i)} - {\overset{\sim}{m}}_{t}^{(j)}}$$M_{t}^{({i,j})}\overset{\Delta}{=}{\min^{2}\left( {{{{d{\overset{\sim}{m}}_{t}^{({i,j})}} - \left( {x + d^{({i,j})}} \right)}},{{{d{\overset{\sim}{m}}_{t}^{({i,j})}} - \left( {x - d^{({i,j})}} \right)}}} \right)}$${\Psi(x)}\overset{\Delta}{=}{\exp\left\{ {- \frac{M_{t}^{({i,j})}}{2\left( \sigma_{ij}^{2} \right)}} \right\}}$$L_{d{\overset{\sim}{m}}_{t}}\overset{\Delta}{=}\left\{ {{\begin{matrix}1 & {{{if}\mspace{14mu}{{x - {d{\overset{\sim}{m}}_{t}^{({i,j})}}}}} < d^{({i,j})}} \\{\Psi(x)} & {otherwise}\end{matrix}{f_{{\overset{\sim}{\Delta}}_{ij}}(x)}}\overset{\Delta}{=}\frac{\prod\limits_{t \in T_{ij}}\; L_{d{\overset{\sim}{m}}_{t}}}{\int_{- \infty}^{\infty}{\prod\limits_{t \in T_{ij}}\; L_{d{\overset{\sim}{m}}_{t}}}}} \right.$where T_(ij) is the set of sources that were measured by both sensor iand sensor j (taken to include all sources, as noted earlier). Thepseudo likelihood function also gives E[{tilde over (Δ)}_(ij)],var[{tilde over (Δ)}_(ij)], and {tilde over (Δ)}_(ij)ϵ(l_(ij), h_(ij))w.p. 1−α.

Collecting the likelihood functions ƒ_({tilde over (Δ)}) _(ij) for alli, j, the {tilde over (δ)} vector can be estimated using the errorcovariance matrix approximation G_({tilde over (δ)})G_({tilde over (δ)})^(T) by using MVIE centering:

Program 3.2 Maximum Volume Inscribed Ellipsoid Centering on ClockOffsets Under Pseudo Likelihood Imposed Bounds

Minimize: log det G_({tilde over (δ)}) ⁻¹

Subject to: ∥G_({tilde over (δ)})a_(i) ^(T)∥+a_(i){tilde over(δ)}_(c)≤b_(i)

-   -   i=1, . . . , M(M−1)

Where:A=[{tilde over (D)} ^(T) −{tilde over (D)} ^(T)]^(T)b=[h _({tilde over (Δ)}) ^(T)(α) l _({tilde over (Δ)})^(T)(α)]^(T)  (3.41)h_({tilde over (Δ)}) ^(T)(α) and l_({tilde over (Δ)}) ^(T)(α) are theupper- and lower-bound vectors imposed by α, via the pseudo likelihoodfunction, on all {tilde over (δ)}pairs.

From {tilde over (δ)}_(c), an estimate of δ can be computed bymultiplying it element-wise by {circumflex over (θ)}:{circumflex over (δ)}={tilde over (δ)}_(c)⊙{circumflex over (θ)}  (3.42)This multiplication cancels the prior normalization and thereforeintroduces no additional errors.

Having estimated {circumflex over (θ)} and {circumflex over (δ)} withtheir respective approximations of the error covariance matrices,G_(θ)G_(θ) ^(T), and G_(δ)G_(δ) ^(T), the skew and bias can be removedfrom the sensor measurements as follows:

$\begin{matrix}\begin{matrix}{{\overset{\sim}{\overset{\sim}{m}}}_{t}^{(s)} = \frac{m_{t}^{s} - {\hat{\delta}}_{s}}{{\hat{\theta}}_{s}}} \\{= {{\frac{\theta_{s}}{{\hat{\theta}}_{s}}\left( {r_{t}^{(s)} + \tau_{t}} \right)} + \frac{\delta_{s} - {\hat{\delta}}_{s}}{{\hat{\theta}}_{s}} + e_{t}^{(s)}}} \\{= {r_{t}^{(s)} + \tau_{t} + {e\;\delta_{s}} - {\tau_{t}e\;\theta_{s}} + e_{t}^{(s)}}}\end{matrix} & (3.43)\end{matrix}$Defining {tilde over ({tilde over (m)})}_(t)=[{tilde over ({tilde over(m)})}_(t) ⁽⁰⁾, {tilde over ({tilde over (m)})}_(t) ⁽¹⁾, . . . , {tildeover ({tilde over (m)})}_(t) ^((M-1))] and calculating the covariance,

$\begin{matrix}{{\Sigma_{t} = {{{var}\left\lbrack {B{\overset{\sim}{\overset{\sim}{m}}}_{t}} \right\rbrack} = {{B\left( {\begin{bmatrix}0 & 0 \\0 & {G_{\delta}G_{\delta}^{T}}\end{bmatrix} + \begin{bmatrix}0 & 0 \\0 & {G_{\theta}G_{\theta}^{T}\tau_{t}^{2}}\end{bmatrix} + {I\;\sigma_{t}^{2}}} \right)}B}}},} & (3.44)\end{matrix}$the optimal source position in the WLS sense reduces to the followingoptimization problem:

$\begin{matrix}{{{argmin}_{p_{t}}\left( {B\left( {{\overset{\sim}{\overset{\sim}{m}}}_{t} - {\frac{1}{c}{{p_{t} - q_{s}}}}} \right)} \right)}^{T}{\Sigma_{t}^{- 1}\left( {B\left( {{\overset{\sim}{\overset{\sim}{m}}}_{t} - {\frac{1}{c}{{p_{t} - q_{s}}}}} \right)} \right)}} & (3.45)\end{matrix}$

Thus, controller 28 can synchronize a sensor network with skewed,asynchronous clocks by measuring the TOA of a set of source signals.Using MVIE to get the initial guess helps avoid local minima in the costfunction.

If the position of one source is known, the clock offset can be measuredat the time of transmission at raw measurement accuracy. If a sourcetransmits twice without changing position, then the clock skew of thesensors that sense it can be estimated directly at the accuracy of theraw measurement. Incorporating such measurements into the system modelcan increase the system accuracy. A source with known location can alsobe a sensor that is allowed to transmit. Periodic sources (such as sometransmitting satellites) usually have a transmission period that is veryshort, so that adjacent transmissions have essentially the sameposition.

As noted earlier, the MVIE estimate together with its accuracyestimation can also be used as a priori probability of the unknowns in aMAP estimator.

The MVIE of the normalized clock offset can be estimated at a few pointsin time to increase accuracy if needed.

Even without assuming that all sources are received by all sensors, theabove algorithm will still perform well with only slight modifications.

Methods Using Combinations of Satellites and Terrestrial Sources

In the description that follows, system 20 is assumed (with fullgenerality) to comprise a set of M sensors 22 at given locations{q_(i)}_(i=0) ^(M-1), a set of N_(s) satellites 26 in orbit of unknownposition {p_(i) ^((s))}_(i=0) ^(N) ^(s) ⁻¹, and a set of N_(t)terrestrial sources 24 of unknown position {p_(i) ^((t))}_(i=0) ^(N)^(t) ⁻¹. Denote by N the total number of sources received by the sensorsN

N_(t)+N_(s), and set the origin of the coordinate system at the centerof the sensor network deployment area. The orientation of the coordinatesystem is set to North West Up.

The time offset of sensor s is δ^((s)). Therefore, the TOA measurementof source t (satellite or terrestrial in origin) by sensor s ism _(t) ^((s)) =r _(t) ^((s))+τ_(t)+δ^((s)) +e _(t) ^((s)),  (4.1)wherein τ_(t) is the unknown transmission time of the signal, r_(t)^((s)) is the time of flight from source t to sensor s, and e_(t) ^((s))is zero mean Gaussian error whose standard deviation is σ_(t) ^((s)).

By defining the vectors:m _(t)

[m _(t) ⁽¹⁾ , m _(t) ⁽²⁾ . . . m _(t) ^((M))]^(T);r _(t)

[r _(t) ⁽¹⁾ , r _(t) ⁽²⁾ . . . r _(t) ^((M))]^(T; r) _(t) ^((s))

1/c∥q _(s) −p _(t)∥1_(M)

[1, 1, . . . 1]^(T); δ

[δ₁, δ₂ . . . δ_(M)]^(T);; e _(t)

[e _(t) ⁽¹⁾ , e _(t) ⁽²⁾ . . . e _(t) ^((M))]^(T).the measurements of source t can be written in vector form as:m _(t) =r _(t)+1_(M)τ_(t) +δ+e _(t).  (4.2)We use the same notation for both terrestrial and satellite sources,wherein the indexes t=0 . . . N_(t)−1 are used for terrestrial sourcesand the indexes t=N_(t) . . . N_(t)+N_(s)−1 are used for satellitesources.

For the sake of simplicity in the formulation we make the assumptions:

1. All of the terrestrial sources are intercepted by all of the sensorswith the same error standard deviation, σ_(T),

2. All of the satellites are intercepted by all of the sensors with thesame error standard deviation σ_(S).

The problem to be solved is twofold:

1. Estimate the sensors clock offsets.

2. Estimate the source location using the available measurements.

For satellite signals the range r_(t) ^((s)) can be written as:

$\begin{matrix}\begin{matrix}{r_{t}^{(s)} = {\frac{1}{c}{{p_{t} - q_{s}}}}} \\{= {\frac{1}{c}\sqrt{{p_{t}}^{2} + {q_{s}}^{2} - {2\; q_{s}^{T}p_{t}}}}} \\{= {\frac{1}{c}\sqrt{{p_{t}}^{2} + {q_{s}}^{2}}\sqrt{1 - \frac{2\; q_{s}^{T}p_{t}}{{p_{t}}^{2} + {q_{s}}^{2}}}}} \\{\approx {\frac{1}{c}{p_{t}}\left( {1 - \frac{q_{s}^{T}p_{t}}{{p_{t}}^{2}} - \frac{\left( {q_{s}^{T}p_{t}} \right)^{2}}{2{p_{t}}^{4}}} \right)}} \\{\approx {\frac{1}{c}{p_{t}}\left( {1 - \frac{q_{s}^{T}p_{t}}{{p_{t}}^{2}}} \right)}} \\{= {{\frac{1}{c}{p_{t}}} - {\frac{1}{c}q_{s}^{T}{\overset{\_}{p}}_{t}}}}\end{matrix} & (4.3)\end{matrix}$

Therefore the vector r_(t) is approximated by

$\begin{matrix}{{r_{t} \approx {{\frac{1}{c}{p_{t}}1_{M}} - {G{\overset{\_}{p}}_{t}}}}{wherein}} & (4.4) \\{G\overset{\Delta}{=}{\frac{1}{c}\left\lbrack {q_{0}\; q_{1}\mspace{14mu}\ldots\mspace{14mu} q_{M - 1}} \right\rbrack}^{T}} & (4.5)\end{matrix}$and p _(t) is a unit vector in the direction of p_(t). The presentmethods use a projection of the range vector that is defined as:

$\begin{matrix}{{{\overset{\sim}{r}}_{t}\overset{\Delta}{=}{Br}_{t}}{B\overset{\Delta}{=}{I - {\frac{1}{M}{11^{T}.}}}}} & (4.6)\end{matrix}$This means that:

$\begin{matrix}{{{\overset{\sim}{r}}_{t} \approx {{B\frac{1}{c}{p_{t}}1_{M}} - {{BG}{\overset{\_}{p}}_{t}}}} = {{- {BG}}{{\overset{\_}{p}}_{t}.}}} & (4.7)\end{matrix}$G is referred to as the geometry matrix. This approximation of theprojected range vector is the reason that the satellites location can bedescribed by two variables. For example, longitude and latitude, insteadof three that are needed for a general point in space.

The second-order term in the Taylor expansion determines theapproximation error. As an example, for GPS satellites at an orbitalheight of 20180 km, and a sensor network of 1 km diameter, the errorterm maximum value is:

$\begin{matrix}{{{\frac{\left( {q_{s}^{T}p_{t}} \right)^{2}}{2{p_{t}}^{3}} \leq \frac{{q_{s}^{T}}^{2}}{2{p_{t}}}} = {\frac{500^{2}}{2 \cdot 20180000} = {6.2\mspace{14mu}{mm}}}},} & (4.8)\end{matrix}$In terms of clock synchronization, this is equivalent to an error of 18ps. Looking at the approximation from the opposite direction, we may askwhat is the maximal radius of the sensor network such that theapproximation error is less then 1 m, for the same satellite orbitalradius. We calculate:

$\begin{matrix}{{\frac{{q_{s}^{T}}^{2}}{2{p_{t}}} = {1\mspace{14mu} m}}{{q_{s}} = {\sqrt{2 \cdot 20180000} = {6352\mspace{14mu}{m.}}}}} & (4.9)\end{matrix}$

Neglecting the second order term is possible, as long as the possibledeployment area of the network is limited. The results can be extendedto wider-area networks, given crude information on the orbital radius ofthe satellites.

Returning to eq. (4.7), the domain of the vector {tilde over (r)}_(t) ispart of a subspace of

^(M). Due to the structure of G, {tilde over (r)}_(t) resides in a two-or three-dimensional subspace of

^(M). The subspace will be two-dimensional for planar sensorconfiguration, and three-dimensional in the general case.

For ground-based sensor configurations, the domain of p _(t) is thethree-dimensional unit half-sphere such that:p _(t) ·z≥0,  (4.10)wherein z is the unit vector pointing up. This inequality holds sincefrom any point on the surface of the earth, the intercepted signalsarrive from above the horizon. In the rare case of space sensorformations, we would have a whole unit sphere as the domain of p _(t).On the other hand, if we have any information on the whereabouts ofsource t, we can use it to reduce the domain of p _(t) to a specifiedsolid angle segment. Although the present description focused onground-based sensor configurations, it can be extended to more generalcases with only slight modifications.

The vector {tilde over (r)}_(t) is a linear transformation of p _(t),and therefore its domain is a linear mapping of the upper half of theunit sphere in

³ onto

^(M), which is a half ellipsoid in

^(M). Singular Value Decomposition (SVD) of −BG reveals the exact shapeof that half ellipsoid. We define the SVD of −BG as:

$\begin{matrix}{{{- {BG}}\overset{\Delta}{=}{{{U\begin{bmatrix}\Lambda_{R} \\0\end{bmatrix}}V^{T}}\overset{\Delta}{=}{{\left\lbrack {R^{T}H^{T}} \right\rbrack\begin{bmatrix}\Lambda_{R} \\0\end{bmatrix}}V^{T}}}},} & (4.11)\end{matrix}$wherein U and V are unitary square matrixes of size M×M, and 3×3respectively; Λ_(R) is the non negative singular value diagonal matrixof size 3×3 at most; R rows hold the left singular vectors related tothe nonzero singular values; and the rest of the left singular vectorsform the rows of H. The nonzero singular values define the lengths ofthe semi-axis of the ellipsoid, while the related left-singular vectorsdefine the semi-axis of the ellipsoid.

When the sensor formation is planar, there are only two non-zerosingular values in Λ_(R). In this case the domain of {tilde over(r)}_(t) is a two-dimensional full ellipsoidal disk in

^(M)

As stated earlier for three-dimensional ground sensor formations, thedomain of p _(t) is a half ellipsoid. To know which half of theellipsoid is the domain, we need to find the plane that separate the twohalves. The plane normal pointing into the domain direction n is:n=−BGz,  (4.12)This plane is not guaranteed to be aligned with the semi-axis of theellipsoid. When the height variation of the sensor configuration is muchsmaller than the planar position variations, however, the plane normalis approximately aligned with the semi-axis of the ellipsoid.

In the above description, there was an implicit condition on the numberof sensors. If there are only two sensors, then rank (BG)=1, and thedomain of BGp _(t) is a line segment. Three sensors are always on aplane or a line, so that the domain BGp _(t) will be a line or anellipsoidal disk. Only four sensors or more can be positioned in anon-planar configuration, and therefore make the domain of BGp _(t) ahalf ellipsoid as discussed.

ML Estimator

The ML estimator for terrestrial sources without clock skew was derivedabove for the localization of a general source using TOA measurements.The cost function to minimize in that case was:

$\begin{matrix}{Q\overset{\Delta}{=}{\sum\limits_{t = 0}^{N - 1}\;{{{m_{t} - r_{t} - {1_{M}\tau_{t}} - \delta}}^{2}.}}} & (4.13)\end{matrix}$where r_(t) is the vector of distances of the source t to all sensors,and τ_(t) is the unknown source transmission time. The configurationcontaining satellites and terrestrial sources of unknown position andtransmission time fits this cost function. Assuming the satellite andterrestrial sources are measured with the same error statistics yields:

$\begin{matrix}{{Q^{\prime}{\sum\limits_{t = 0}^{N_{T} - 1}{{m_{t} - {r_{t}\left( {x_{t},y_{t},z_{t}} \right)} - {1_{M}\tau_{t}} - \delta}}^{2}}} + {\sum\limits_{t = N_{T}}^{N_{T} + N_{S} - 1}{{m_{t} - {r_{t}\left( {\phi_{t},\theta_{t}} \right)} - {1_{M}\tau_{t}} - \delta}}^{2}}} & (4.14)\end{matrix}$where r_(t) (x_(t), y_(t), z_(t)) is the distance vector for aterrestrial source located at p_(t)

[x_(t), y_(t), z_(t)]^(T), and r_(t) (φ_(t), θ_(t)) is the distancevector for a satellite source at azimuth θ_(t) and elevation φ_(t). Thetransmission times τ_(t) that minimize the cost function are given by:

$\begin{matrix}\begin{matrix}{{\hat{\tau}}_{t} = {\left( {1_{M}^{T}1_{M}} \right)^{- 1}1_{M}^{T}\left( {m_{t} - r_{t} - \delta} \right)}} \\{= {\frac{1}{M}1_{M}^{T}{\left( {m_{t} - r_{t} - \delta} \right).}}}\end{matrix} & (4.15)\end{matrix}$

We now define:

$\begin{matrix}{{{{BI_{M}} - {\frac{1}{M}1_{M}1_{M}^{T}}};}{{{\overset{\sim}{m}}_{t}{Bm}_{t}};{{\overset{\sim}{r}}_{t}{Br}_{t}};{\overset{\sim}{\delta}B\;{\delta.}}}} & (4.16)\end{matrix}$Substituting the transmission time estimates into the cost function andusing these definitions gives:

$\begin{matrix}{Q^{\prime} = {{\sum\limits_{t = 0}^{N_{T} - 1}{{{\overset{\sim}{m}}_{t} - {{\overset{\sim}{r}}_{t}\left( {x_{t},y_{t},z_{t}} \right)} - \overset{\sim}{\delta}}}^{2}} + {\sum\limits_{t = N_{T}}^{N_{T} + N_{s} - 1}{{{\overset{\sim}{m}}_{t} - {r_{t}\left( {\phi_{t},\theta_{t}} \right)} - \overset{\sim}{\delta}}}^{2}}}} & (4.17)\end{matrix}$Collecting all vectors for t=0, 1, . . . , N_(S)+N_(T)−1 now gives:Q′=∥m−r−Cδ∥ ²,  (4.18)wherem=

[{tilde over (m)} ₀ ^(T) , {tilde over (m)} ₁ ^(T) . . . {tilde over(m)} _(N) _(T) _(+N) _(s) ⁻¹ ^(T) ; r

[{tilde over (r)} ₀ ^(T) , {tilde over (r)} ₁ ^(T) . . . {tilde over(r)} _(N) _(T) _(+N) _(s) ⁻¹]^(T);e

[{tilde over (e)} ₀ ^(T) , {tilde over (e)} ₁ ^(T) . . . {tilde over(e)} _(N) _(T) _(+N) _(s) ⁻¹]^(T) ; C

[I _(M) , I _(M) . . . I _(M)]^(T).  (4.19)

Substituting the clock offsets into the cost function gives:

$\begin{matrix}{Q^{\prime} = {{{\left\lbrack {I - {\frac{1}{N}{CC}^{T}}} \right\rbrack\left( {m - r} \right)}}^{2}.}} & (4.20)\end{matrix}$Calculating the elements of r is straightforward for terrestrialsources. For satellite sources, we can use the fact that they are faraway to write:

$\begin{matrix}{{{{\overset{\sim}{r}}_{t}{\lim\limits_{R->\infty}{Br}_{t}}} = {{- {BG}}{\overset{\_}{p}}_{t}}},} & (4.21)\end{matrix}$where R is the orbital radius of the satellite sources. Thus, findingthe ML estimate requires a search in 3N_(T)+2N_(S) dimensions: threeposition coordinates for each terrestrial source and two bearing anglesfor each celestial source. As noted above, for ground-based sensorformations, the satellite directions are above the horizon. Therefore:

$\begin{matrix}{\phi_{t} \in {\left( {0,\frac{\pi}{2}} \right).}} & (4.22)\end{matrix}$This criterion limits the search interval.

To avoid a brute-force search over so many dimensions, an estimate of{tilde over (δ)} can be used to decouple the problems of estimating theterrestrial source positions and satellite directions. The followingiterative procedure can then be used to solve the ML estimator:

1. Initialize: {tilde over (δ)}⁽⁰⁾←0; Q′⁽⁻¹⁾←inf; n←0

2. Q′^((n))←0

3. For each source t=1, 2, . . . , N:

-   -   (a) Estimate p_(t) given {tilde over (δ)}^((n)) by minimizing        ∥{tilde over (m)}_(t)−{tilde over (r)}_(t)−{tilde over        (δ)}^((n))∥²    -   (b) Q′^((n))←Q′^((n))+∥{tilde over (m)}_(t)−{tilde over        (r)}_(t)({circumflex over (p)}_(t))−{tilde over (δ)}^((n))∥²

4.

$\left. {\overset{\sim}{\delta}}^{({n + 1})}\leftarrow{\frac{1}{N}{{C^{T}\left( {m - r} \right)}.}} \right.$

5. If Q′_(Δ) ^((n))=|Q′^((n))−Q′^((n−1))|<ϵ then stop

6. Else n←n+1, go to step 2.

If the computation is found to have converged (step 56 in FIG. 4),controller 28 synchronizes the clocks of sensors 22 based on the currentvalue of {tilde over (δ)} (step 58). Depending on applicationrequirements, the controller may either send actual synchronizationvalues to the sensors, or it may simply record and apply appropriatecorrections to the clock readings received from the sensors. Thecontroller may also output or otherwise apply the source location valuesgiven by the current elements of r (step 60).

For every iteration of this algorithm, controller 28 performsN_(S)+N_(T) nonlinear least square optimizations. To avoid convergenceof the cost function to a local minimum (rather than the true minimumcorresponding to the actual clock offsets and source positions), a goodinitial estimate of the source locations and/or clock offsets is useful.Such an estimate can reduce the computational load considerably byreducing the number of iterations required for convergence. In trackingapplications, a previous estimate may be used as a good starting point,but in the absence of a previous estimate, other information may beused.

The following sections describe methods that can be used to find aninitial point so that no more than a few iterations are needed.Moreover, given a good initial point, controller 28 can solve theoriginal cost function using gradient methods, which are very effective.

Locating Satellite Direction

In the iterative method of FIG. 4, the direction of a satellite given aguess of the clock offset vector is required. In this section we presenta Lagrange multiplier approach to finding the required direction, whichmay be preferable to grid search methods. For this purpose, theoptimization of a set of satellite source locations can be written asfollows:

$\begin{matrix}{Minimize} & {{{Bm}_{t} + {{BG}{\overset{\_}{p}}_{t}} - {B\;\delta}}}^{2} \\{\overset{\_}{p}}_{t} & \; \\{{Subject}\mspace{14mu}{to}} & {{{\overset{\_}{p}}_{t}} = 1.}\end{matrix}$We define:y

Bδ−Bm _(t)A

BG.  (4.24)Then the Lagrangian is:

$\begin{matrix}{{L = {{{{A{\overset{\_}{p}}_{t}} - y}}^{2} - {\lambda\left( {1 - {{\overset{\_}{p}}_{t}^{T}{\overset{\_}{p}}_{t}}} \right)}}}{{\frac{\partial}{\partial{\overset{\_}{p}}_{t}}L} = {{A^{T}\left( {{A{\overset{\_}{p}}_{t}} - y} \right)} + {\lambda{\overset{\_}{p}}_{t}}}}{{\frac{\partial}{\partial\lambda}L} = {1 - {{{\overset{\_}{p}}_{t}}^{2}.}}}} & (4.25)\end{matrix}$

Equating

$\frac{\partial}{\partial{\overset{\_}{p}}_{t}}$L to zero gives:(A ^(T) A−Iλ) p _(t) =A ^(T) yp _(t)(A ^(T) A−Iλ)⁻¹ A ^(T) y.   (4.26)Substituting the solution for p _(t) in the expression for

$\frac{\partial}{\partial\lambda}L$and equating to zero gives:y ^(T) A(A ^(T) A−Iλ)⁻¹ ^(T) (A ^(T) A−Iλ)⁻¹ A ^(T) y=1.   (4.27)A^(T)A can be represented by SVD as:A ^(T) A=V ^(T) ΛVy ^(T) AV ^(T)(Λ−Iλ)⁻¹(Λ−Iλ)⁻¹ VA ^(T) y=1.  (4.28)

Noting that Λ

diag (γ₀, γ₁, γ₂), and defining VA^(T)y

α=[α₀, α₁, α₂]^(T) we get:

$\begin{matrix}{{\frac{\alpha_{0}^{2}}{\left( {\gamma_{0} - \lambda} \right)^{2}} + \frac{\alpha_{1}^{2}}{\left( {\gamma_{1} - \lambda} \right)^{2}} + \frac{\alpha_{2}^{2}}{\left( {\gamma_{2} - \lambda} \right)^{2}}} = 1.} & (4.29)\end{matrix}$This equation is equivalent to finding the roots of a sixth degreepolynomial. An efficient numerical solution for the roots can be foundby seeking the eigenvalues of the companion matrix of the polynomial.One way to reduce the complexity in runtime is by realizing that the γvalues are constant as long as the sensors are static, they are neededonly once, at initialization.

After calculating the six possible values of λ_(opt) and discarding anycomplex solutions, we place them in the equation for p _(t) to getpossible directions of the satellite. Controller 28 can test which ofthe real solutions is the optimal one by evaluating the cost functionfor each possible value, thus assuring that the global optimum ischosen.

A weight matrix W may be used in the original problem:

$\begin{matrix}\begin{matrix}{Minimize} & {{W\left( {{Bm}_{t} + {{BG}{\overset{\_}{p}}_{t}} - {B\;\delta}} \right)}}^{2} \\{\overset{\_}{p}}_{t} & \; \\{Subjectto} & {{{{\overset{\_}{p}}_{t}} = 1},}\end{matrix} & (4.30)\end{matrix}$To solve the problem, y is replaced with Wy and A with WA. The use ofthe weight matrix extends the results to include correlated measurementsand can be applied to reduce the dimension of the measurement vector, asexplained below.

ML Decomposition for Satellite-Transmitted Signals

This section shows that the log likelihood cost function for signalsfrom space can be divided into two separate cost functions. One islinear, while the other is nonlinear but is only three-dimensional.Controller 28 may decompose the cost function into such components inorder to solve the problems of clock synchronization and satellitesource direction more efficiently.

H and R are defined above as the null and range-space projectors of BGp_(t). We define:T

[H ^(T) , R ^(T)]^(T).  (4.31)The matrix T is a unitary matrix, since it is a square matrix withorthonormal rows. Multiplying the source measurements m_(t) by T doesnot reduce the information in the measurements. Proceeding with themeasurement of the satellite signals gives:Tm _(t) =Tr _(t) +T1_(M)τ_(t) +Tδ+Te _(t)  (4.32)

Solving the ML estimator then requires that controller 28 minimize thefollowing cost function:

$\begin{matrix}{Q{\sum\limits_{t = 0}^{N_{s} - 1}{{{T\left( {m_{t} - r_{t} - {1_{M}\tau_{t}} - \delta} \right)}}^{2}.}}} & (4.33)\end{matrix}$Transmission time is estimated as usual resulting in:

$\begin{matrix}\begin{matrix}{Q = {\sum\limits_{t = 0}^{N_{s} - 1}\;{{\left\lbrack {H^{T},R^{T}} \right\rbrack^{T}\left( {{Bm}_{t} - {Br}_{t} - {B\;\delta}} \right)}}^{2}}} \\{= {{\sum\limits_{t = 0}^{N_{s} - 1}\;{{H\left( {{Bm}_{t} - {B\;\delta}} \right)}}^{2}} + {\sum\limits_{t = 0}^{N_{s} - 1}\;{{R\left( {{Bm}_{t} - {Br}_{t} - {B\;\delta}} \right)}}^{2}}}}\end{matrix} & (4.34)\end{matrix}$since HBr_(t)=0 by definition of H.

We now define:δ^(⊥)

HBδδ^(R)

RBδ.  (4.35)Therefore:

$\begin{matrix}{Q = {{\sum\limits_{t = 0}^{N_{s} - 1}\;{{{HBm}_{t} - \delta^{\bot}}}^{2}} + {\sum\limits_{t = 0}^{N_{s} - 1}\;{{{{RBm}_{t} - {RBr}_{t} - \delta^{R}}}^{2}.}}}} & (4.36)\end{matrix}$

The two parts of the cost Q can be solved independently, thereby makingthe ML estimator simpler to compute. Because the measurement vectorsHBm_(t) and RBm_(t) are uncorrelated, the estimation error covariancematrix is block diagonal. The ML estimators related to the two costfunctions can thus be computed separately, as described below.

The ML estimator of δ^(⊥) can be derived from the partial cost function:

$\begin{matrix}{Q^{\bot}\overset{\Delta}{=}{\sum\limits_{t = 0}^{N_{s} - 1}\;{{{{HBm}_{t} - \delta^{\bot}}}^{2}.}}} & (4.37)\end{matrix}$For this purpose, we define:m _(t) ^(⊥)

HBm _(t)e _(t) ^(⊥)

HBe _(t)m ^(⊥)

[m ₀ ^(⊥) ^(T) , m ₁ ^(⊥) ^(T) , . . . , m _(N) _(s) ⁻¹ ^(⊥) ^(T) ]^(T)e ^(⊥)

[e ₀ ^(⊥) ^(T) , e ₁ ^(⊥) ^(T) , . . . , e _(N) _(s) ⁻¹ ^(⊥) ^(T) ]^(T)C ^(⊥)

1_(N) _(s)

I  (4.38)wherein

is the Kronecker Product. Now we can write:Q ^(⊥) =∥m ^(⊥) −C ^(⊥)δ^(⊥)∥².  (4.39)The estimator of δ^(⊥) that minimizes the cost function is:

$\begin{matrix}{{\hat{\delta}}^{\bot} = {\frac{1}{N}C^{\bot^{T}}{m^{\bot}.}}} & (4.40)\end{matrix}$

It can be shown that HBH^(T) is idempotent, and H can be chosen suchthat:

$\begin{matrix}{{HBH}^{T} = {\begin{bmatrix}0 & 0 \\0 & I\end{bmatrix}.}} & (4.41)\end{matrix}$

Furthermore, regardless of the source locations, the estimation accuracydepends only on the number of signals measured. The positions of sensors22 also do not affect the estimation accuracy, as long as the sensorsare confined to an area such that the linear approximation of thesatellite ranges holds.

A single source is sufficient to establish an estimation of the nullspace projection of the clock offsets (although this is not the case forthe range space). Estimating the clock offsets in null space issufficient to localize terrestrial sources with reasonable accuracy. Fora non-planar sensor configuration, sources confined to a plane can belocalized if at least six sensors observe the signals, and for sourcesthat are not in a plane at least seven sensors are required. For aplanar sensor configuration, sources confined to a plane can belocalized if at least five sensors observe the signals, whereas forsources that are not confined to a plane, at least six sensors arerequired.

The fact that we can choose H such that HBH^(T) has zero off-diagonalelements can speed up the estimation process.

For satellite signals, the nonlinear part of the cost function can beexpressed as follows:

$\begin{matrix}{Q^{R} = {\sum\limits_{t = 0}^{N_{s} - 1}\;{{{{RBm}_{t} - {RBr}_{t} - \delta^{R}}}^{2}.}}} & (4.42)\end{matrix}$Controller 28 can estimate δ^(R) in the presence of RBr_(t) as anuisance parameter. In this calculation, since R=RB, the measurementvector RBm_(t) elements are independent and identically distributedGaussian random variables. The mode of calculation depends on the sensorformation: If the sensors are confined to a plane, then RBr_(t) can beany value inside an ellipse, whereas if the sensor distribution is in 3Dspace, RBr_(t) is distributed on a 3D ellipsoidal half shell, whichenables better accuracy.

First we consider the case of a planar sensor configuration. In thiscase the problem of estimating δ^(R) is reduced to finding the possiblecenters of an ellipse of known dimensions given a set of points in theellipse domain corrupted by noise. The set of possible centers can befound, for example, by calculating the mean of the points (which isguaranteed to be inside of the ellipse) and then expanding the set ofpossible points by propagating out from that point. Controller 28 candetermine whether a point is in the domain of possible offsets bychecking whether an elliptical disk with its center at that pointinclude all of the measured points in its domain. Regardless of how theset of points is found, its mean can be written as {circumflex over(δ)}^(R), and its empirical covariance matrix as Σ_(R).

For the case of a non-planar sensor configuration, because RBr_(t)resides on a 3D ellipsoidal half-shell, estimating {circumflex over(δ)}^(R) is reduced to finding the ellipsoid center given a set ofpoints that are located on the ellipsoid half shell in the noiselesscase, or are near this shell in the noisy case. This problem can besolved in a number of ways, including high-dimensional ML estimation,low-dimensional iterative ML estimation, and ellipsoid surface fitting.

The high-dimensional ML estimator takes the range space cost functionand minimizes it for all satellite directions p _(t) and range-projectedclock offsets δ^(R). This is an optimization problem in 3+2N_(S)dimensions which is easy to solve as long as the number of satellites issmall enough.

The low-dimensional iterative ML estimator is similar to the iterativeML algorithm described above, but is applied only to minimize the costfunction with respect to satellite directions and a 3D δ^(R) vector.Each iteration solves two or three dimensions at a time and isinherently fast. This method uses a two-step iterative process: Untilconvergence do ML Estimation of δ^(R) given the current estimation of p_(t); and then update the values of {p _(t)}_(t=0) ^(N) ^(s) ⁻¹ using MLestimation based on the new δ^(R) estimate. This step can beaccomplished by the Lagrange multiplier technique presented above, usingR as the weighting matrix W. This estimation process thus comprises anaveraging operation followed by an eigenvalue calculation of a 6×6matrix. This estimator is therefore fast as long as the iteration countis reasonable.

Estimation Clock Offset Range Space Using Ellipsoid Fitting

The surface fitting method uses the fact that an ellipsoid center is tobe fitted to noisy points, which is again a 3D problem with fastconvergence properties. It is based generally on techniques described byAhn et al., in “Orthogonal Distance Fitting of Implicit Curves andSurfaces,” IEEE Transactions on Pattern Analysis and MachineIntelligence, 24(5):620-638 (2002), which is incorporated herein byreference. This fitting method is a good approximation of the MLestimator and is relatively simple to apply in the present case sincethe ellipsoid to be fitted is of known shape and orientation, so theonly parameters left to fit are the ellipsoid center coordinates.

The problem of 3D ellipsoid fitting can be framed as follows: Given theset of noisy points {RBm_(t)}_(t=0) ^(N) ^(s) ⁻¹, find the 3D ellipsoidcenter that minimizes the weighted sum of the squared distances from thenoisy measurements to the ellipsoid. The cost function to minimize isthus:C(a)=d ^(T) P ^(T) Pdd=[d ₁ , d ₂ , . . . d _(N) _(s) ]d _(i) =∥x _(i) −x′ _(i)(a)∥  (4.43)Here a is the proposed ellipsoid center, x is the noisy point, x′ is theclosest point to x that is on the ellipsoid centered at a, P^(T)P is aweighting matrix, and d_(i) is the Euclidean distance between x_(i) andx′_(i). It is easier to consider x, x′ in a coordinate system that iscentered at a and aligned with the eigenvectors given by the rows of R.In that coordinate system the ellipsoid is given by the followingimplicit function:

$\begin{matrix}{{F\left( x^{\prime} \right)} = {{\frac{x_{1}^{\prime 2}}{\lambda_{1}} + \frac{x_{2}^{\prime 2}}{\lambda_{2}} + \frac{x_{3}^{\prime 2}}{\lambda_{3}} - 1} = 0}} & (4.44)\end{matrix}$wherein x′=[x′₁, x′₂, x′₃]^(T), and λ_(i) i=1, 2,3 are the non-zeroeigenvalues of BGG^(T)B^(T), related to the eigenvectors defined by therows of R.

The ellipsoid center estimation may be carried out using the followingorthogonal distance fitting algorithm to find the ellipsoid center â:

1. â←â₀ Δa←∞k←0

2. While ∥Δa∥>v_(thr)

-   -   2.1. For All measurements x, calculate x′ using Eq. (4.45)    -   2.2. Calculate Δa using Eq. (4.47)    -   2.3. â←â+η^(k)Δa    -   2.4. k←k+1

End While

3. Calculate cov(â) using Eq. (4.48)

4. δ^(R)←â

5. Σ^(R)←cov(â)

In step 2.2, x′ can be found using a generalized Newton-Raphson method,starting from the initial guess x′₀=x:

$\begin{matrix}{{\left. \frac{\partial f}{\partial x^{\prime}} \middle| {}_{k}{\Delta\; x^{\prime}} \right. = \left. {- {f(x)}} \right|_{k}},{x_{k + 1}^{\prime} = {x_{k}^{\prime} + {\Delta\; x^{\prime}}}}} & (4.45)\end{matrix}$wherein:

$\begin{matrix}{f = \begin{bmatrix}{\frac{x_{1}^{\prime 2}}{\lambda_{1}} + \frac{x_{2}^{\prime 2}}{\lambda_{2}} + \frac{x_{3}^{\prime 2}}{\lambda_{3}} - 1} \\{{{- 2}\; x_{1}{\lambda_{2}\left( {x_{2} - x_{2}^{\prime}} \right)}} + {2\; x_{2}{\lambda_{1}\left( {x_{1} - x_{1}^{\prime}} \right)}}} \\{{{- 2}\; x_{2}{\lambda_{3}\left( {x_{3} - x_{3}^{\prime}} \right)}} + {2\; x_{3}{\lambda_{2}\left( {x_{3} - x_{3}^{\prime}} \right)}}}\end{bmatrix}} & (4.46)\end{matrix}$

In step 2.2, Δa is found using:J ^(T) P ^(T) PJΔa=−J ^(T) P ^(T) PdΔa=−(J ^(T) P ^(T) PJ)⁻¹ J ^(T) P ^(T) Pd| _(k)â _(k+1) =â _(k)+η^(k) Δa  (4.47)wherein

${J\overset{\Delta}{=}\frac{\partial d}{\partial a}},$and ηϵ(0,1) is a constant that helps convergence in noisy situations. Insimulations, the inventors found it useful to set η=0.9.

In step 3, the covariance is found using{circumflex over (δ)}R=âΣ_(R) =cov({circumflex over (δ)}^(R))=σ_(S) ²(J ^(T) J)⁻¹.  (4.48)

The above algorithm for ellipsoid center estimation requires invertingmatrices only of size 3×3. Because this algorithm uses the steepestdescent at the ellipsoid center update, as defined by Δa, whereas theprevious algorithm uses ML estimation of the center under the assumptionthat the satellites directions are known, the present algorithm mayconverge in fewer steps, although each step is a bit more complex.

Regardless of the method used to estimate the range space, it isimportant to start the iterative process with a good initialization of â(or {circumflex over (δ)}^(R)) in order to avoid reaching a localminimum. The points to be fitted, assuming a small sensor heightvariation, reside on the upper or lower half of an ellipsoid, and thesign of −RBG·z can be evaluated to identify the part of the ellipsoid onwhich the points are located. (If the height variation is significant,then the initial guess is different but can still be derived in astraightforward manner.) A good initial guess for the ellipsoid center,â, would put the measured points in the vicinity of the correct half ofthe ellipsoid. Such an initial guess could be performed as follows:

1. Set (RBm) _(xy) equal the average for each of the x,y coordinates ofRBm_(t).

2. â₀ _(xy) ←(RBm) _(xy)

3. Define: ϵ>0

4. if −RBG·{circumflex over (z)}≥0

-   -   4.1 â₀ _(z) ←min((RBm)_(z))−ϵ

5. Else

-   -   5.1 â₀←max((RBm)_(z))+ϵ

EndIf

6.

$\left. {\hat{a}}_{0}\leftarrow\begin{bmatrix}{\hat{a}}_{0_{xy}}^{T} & {\hat{a}}_{0_{z}}\end{bmatrix}^{T} \right.$

Fusing Null and Range Space Clock Offset Estimates

By combining the above estimates for the null space and range space ofthe clock offset vector Bδ, controller 28 can now establish fullsynchronization of sensors 22 using satellite sources 26. For thispurpose, the estimations of the range and null projections may betreated as the measurements:

$\begin{matrix}{\begin{bmatrix}{\hat{\delta}}^{\bot} \\{\hat{\delta}}^{R}\end{bmatrix} = {{{\begin{bmatrix}H \\R\end{bmatrix}B\;\delta} + \begin{bmatrix}{nR} \\{nH}\end{bmatrix}}\overset{\Delta}{=}{{{TB}\;\delta} + {n.}}}} & (4.49)\end{matrix}$The error elements n_(H), n_(R) have the covariance matricesΣ_({circumflex over (δ)}) _(⊥) and Σ^(R) respectively and areuncorrelated. Since Σ_({circumflex over (δ)}) _(⊥) is singular, the MLsolution is given by:

$\begin{matrix}{= {\left( {T^{T}\Sigma^{\dagger}T} \right)^{\dagger}T^{T}{\Sigma^{\dagger}\begin{bmatrix}{\hat{\delta}}^{\bot} \\{\hat{\delta}}^{R}\end{bmatrix}}}} & (4.50)\end{matrix}$wherein Σ

blkdiag (Σ_({circumflex over (δ)}) _(⊥) , Σ^(R)), and blkdiag is a blockmatrix whose blocks are the indicated matrices. Since T is unitary, theML estimator becomes:

$\begin{matrix}{\begin{matrix}{= {T^{T}\Sigma\;{TT}^{T}{\Sigma^{\dagger}\begin{bmatrix}{\hat{\delta}}^{\bot} \\{\hat{\delta}}^{R}\end{bmatrix}}}} \\{= {T^{T}\begin{bmatrix}{\hat{\delta}}^{\bot} \\{\hat{\delta}}^{R}\end{bmatrix}}}\end{matrix}{{{cov}{()}} = {T^{T}\Sigma\;{T.}}}} & (4.51)\end{matrix}$

To estimate the directions of the satellite sources,

is inserted into the ML cost function:Q _(t)

∥B(m _(t) −r _(t))−

∥²,  (4.52)and controller 28 finds the solution that minimizes the cost for everysatellite t. Alternatively, the Lagrange multipliers described above maybe used for this purpose.

Thus, controller 28 is able to synchronize a set of sensors 22 deployedover a wide area using satellites having unknown orbital parameters andtime references. The directions to the satellites are estimated as aby-product. The accuracy of the synchronization is only limited by theradio bandwidth, signal/noise ratio (SNR), observation time, and thenumber of sources processed.

Locating Terrestrial Sources

Given only the null-space offset projection {circumflex over (δ)}^(⊥),controller 28 can establish synchronized time-based localization ofterrestrial sources 24. This approach can be useful since estimating{circumflex over (δ)}^(R) is more demanding, and cannot be accomplishedby the above methods with fewer than three celestial signals. Moreover,this sort of localization enables estimation of the full clock offsetvector Bδ, as explained below.

Although HBr_(t)=0 for celestial sources, this is not the case forterrestrial sources 24. Therefore, if there is a sufficient number ofsensors 22 receiving terrestrial signals, it is possible to get asynchronized localization of the terrestrial sources by solving:HBm _(t) =HBr _(t) +HBδ+HBe _(t)HBm _(t)−{circumflex over (δ)}^(⊥) =HBr _(t) +e _({circumflex over (δ)})_(⊥) +HBe _(t)  (4.53)wherein m_(t) is now the vector of TOA measurements of terrestrialsource t, and e_({circumflex over (δ)}) _(⊥) is the estimation error of{circumflex over (δ)}^(⊥).

Since the error e_({circumflex over (δ)}) _(⊥) is uncorrelated withe_(t), the expected estimation error under small error approximation isgiven by:

$\begin{matrix}{\begin{matrix}{{{cov}\left( p_{t} \right)} = \left( {L_{t}^{T}H^{T}\Sigma_{t}^{\dagger}{HL}_{t}} \right)^{- 1}} \\{= {\left( {\frac{\sigma_{S}^{2}}{N_{S}} + \sigma_{T}^{2}} \right)\left( {L_{t}^{T}H^{T}{HBH}^{T}{HL}_{t}} \right)^{- 1}}} \\{{= {\left( {\frac{\sigma_{S}^{2}}{N_{S}} + \sigma_{T}^{2}} \right)\left( {{L_{t}^{T}L_{t}} - {L_{t}^{T}R^{T}{RL}_{t}}} \right)^{- 1}}},}\end{matrix}{wherein}{L_{t}\overset{\Delta}{=}{\frac{\partial{Br}_{t}}{\partial p_{t}} = {\frac{\partial{\overset{\sim}{r}}_{t}}{\partial p_{t}}.}}}} & (4.54)\end{matrix}$

Alternatively, controller 28 may locate terrestrial sources 24 using thefully-estimated offset vector and the terrestrial source sensormeasurements:Bm _(t) =Br _(t) +Bδ+Be _(t)  (4.55)In this case, p_(t) is found by minimizing:argmin_(p) _(t) (Bm _(t) −Br _(t)−

)^(T){tilde over (Σ)}_(t) ^(†)(Bm _(t) −Br _(t)−

){tilde over (Σ)}_(t) =T ^(T) ΣT+Bσ _(T) ².  (4.56)The expected estimation error under linear approximation is given by:cov(p _(t))=(L _(t) ^(T){tilde over (Σ)}_(t) ^(†) L _(t))⁻¹  (4.57)This estimation will thus provide better localization performance thanthe solution based only on null-space offset projection.

Refining Clock Offset Estimation Using Terrestrial Source

The accuracy of the clock offset estimation for sensors 22 may beimproved by considering the localization of terrestrial sources 24.After refining the clock offsets, controller 28 can improve thelocalizations of the terrestrial sources and continue iterating toconvergence. Convergence is assured, as improving the clock offsetestimation accuracy increases the accuracy of localization of theterrestrial sources and of the directions of the satellite sources, andvice versa.

The ML estimator for Bδ may be derived in this case as follows: Theerror of {circumflex over (p)}_(t) is a vector in two or threedimensions, whereas {tilde over ({circumflex over (r)})} is a vector inM dimensions. We define the error vectors:e _(p) _(t)

p _(t) −{circumflex over (p)} _(t)e _({tilde over (r)}) _(t)

{tilde over (r)} _(t)−{tilde over ({circumflex over (r)})}_(t).  (4.58)

There is a linear relationship between the error vectors, assuming theerror is small:e _({tilde over (r)}) _(t) =L _(t) e _(p) _(t) .  (4.59)The error covariance matrix of {tilde over (r)}_(t) is given by:cov({tilde over ({circumflex over (r)})}_(t))=L _(t) cov(p _(t))L _(t)^(T),  (4.60)which is clearly singular and of the same rank as cov(p_(t)) (which isrank 2 or 3 depending on the sources distribution in a plane or in threedimensional space). Eigenvalue decomposition of cov({tilde over({circumflex over (r)})}) gives:cov({tilde over ({circumflex over (r)})}_(t))=U _(t)Λ_(t) U _(t) ^(T)U _(t) =[R _(t) ^(T) , H _(t) ^(T)]Λ_(t)=diag(λ_(t) ₁ , λ_(t) ₂ , λ_(t) ₃ , 0, . . . , 0)R _(t) R _(t) ^(T) =IH _(t) H _(t) ^(T) =I,  (4.61)where the rows of R_(t) are the eigenvectors of cov({tilde over({circumflex over (r)})}_(t)) corresponding to the non-zero eigenvalues,while the rows of H_(t) are the eigenvectors corresponding to the zeroeigenvalue, and λ_(t) ₃ =0 in the case of planar sensor formation.

Now we calculate the covariance of H_(t)e_({tilde over (r)}) _(t) :

$\begin{matrix}\begin{matrix}{{{cov}\left( {H_{t}e_{{\overset{\sim}{r}}_{t}}} \right)} = {H_{t}U_{t}\Lambda_{t}U_{t}^{T}H_{t}^{T}}} \\{= {{H_{t}\left\lbrack {R_{t}^{T},H_{t}^{T}} \right\rbrack}{\Lambda_{t}\left\lbrack {R_{t}^{T},H_{t}^{T}} \right\rbrack}^{T}H_{t}^{T}}} \\{= {\left\lbrack {0,I} \right\rbrack{\Lambda_{t}\left\lbrack {0,I} \right\rbrack}}} \\{= 0}\end{matrix} & (4.62)\end{matrix}$This expression means that the first two moments of the projected errorH_(t)e_({tilde over (r)}) _(t) are zero, and under the assumption ofGaussian distribution yields:H _(t) e _({tilde over (r)}) _(t) =0, Almost Sure.  (4.63)Since p_(t) is unknown, the evaluation of L_(t) and cov(p_(t)) isperformed using the estimate of p_(t).

Coming back to the source measurement equation, multiplying from theleft by H_(t) gives:H _(t) Bm _(t) =H _(t) {tilde over (r)} _(t) +H _(t) Bδ+H _(t) Be_(t)  (4.64)This equation is modified to compensate by {tilde over ({circumflex over(r)})}_(t):

$\begin{matrix}\begin{matrix}{{H_{t}\left( {{Bm}_{t} - {\hat{\overset{\sim}{r}}}_{t}} \right)} = {{H_{t}e_{{\overset{\sim}{r}}_{t}}} + {H_{t}B\;\delta} + {H_{t}{Be}_{t}}}} \\{{= {{H_{t}B\;\delta} + {H_{t}{Be}_{t}}}},}\end{matrix} & (4.65)\end{matrix}$The result is a linear measurement of a projection of the sensor clockoffsets vector. Concatenating measurements from all of the terrestrialsources gives:

$\begin{bmatrix}{H_{0}\left( {{Bm}_{0} - {\hat{\overset{\sim}{r}}}_{0}} \right)} \\{H_{1}\left( {{Bm}_{1} - {\hat{\overset{\sim}{r}}}_{1}} \right)} \\\vdots \\{H_{N_{t} - 1}\begin{pmatrix}{{Bm}_{N_{t} - 1} -} \\{\hat{\overset{\sim}{r}}}_{N_{t} - 1}\end{pmatrix}}\end{bmatrix} = {{\begin{bmatrix}H_{0} \\H_{1} \\\vdots \\H_{N_{t} - 1}\end{bmatrix}B\;\delta} + {\begin{bmatrix}{H_{0}B} & \; & \; & 0 \\\; & {H_{1}B} & \; & \; \\\; & \; & \ddots & \; \\0 & \; & \; & {H_{N_{t} - 1}B} \\\; & \; & \; & \;\end{bmatrix}\begin{bmatrix}e_{0} \\e_{1} \\\vdots \\{e_{N_{t} - 1},}\end{bmatrix}}}$which is completely independent of the localization errors.

Writing all of the equations we have on Bδ at this point gives:

$\begin{bmatrix}{H_{0}\left( {{Bm}_{0} - {\hat{\overset{\sim}{r}}}_{0}} \right)} \\{H_{1}\left( {{Bm}_{1} - {\hat{\overset{\sim}{r}}}_{1}} \right)} \\\vdots \\{H_{N_{t} - 1}\begin{pmatrix}{{Bm}_{N_{t} - 1} -} \\{\hat{\overset{\sim}{r}}}_{N_{t} - 1}\end{pmatrix}} \\ \\

\end{bmatrix} = {\quad{\begin{bmatrix}H_{0} \\H_{1} \\\vdots \\H_{N_{t} - 1} \\H \\R\end{bmatrix}{\quad{B\;\delta{\quad{\quad{\quad{+ {\quad{\begin{bmatrix}{H_{0}B} & \; & \; & \; & \; & \; \\\; & {H_{1}B} & \; & \; & 0 & \; \\\; & \; & \ddots & \; & \; & \; \\\; & \; & \; & {H_{N_{t} - 1}B} & \; & \; \\\; & {0\;} & \; & \; & {\sqrt{\frac{\sigma_{S}^{2}}{N_{S}}}{HB}} & \; \\\; & \; & \; & \; & \; & \sqrt{\sum\limits_{R}^{\;}}\end{bmatrix}\begin{bmatrix}e_{0} \\e_{1} \\\vdots \\e_{N_{t}} \\e^{\bot} \\e_{R}\end{bmatrix}}\;}}}}}}}}}$

The ML estimator is established using WLS with the followingdefinitions:

$m_{B\;\delta}\overset{\bigtriangleup}{=}\begin{bmatrix}\left( {H_{0}\left( {{Bm}_{0} - {\hat{\overset{\sim}{r}}}_{0}} \right)} \right)^{T} & \left( {H_{1}\left( {{Bm}_{1} - {\hat{\overset{\sim}{r}}}_{1}} \right)} \right)^{T} & \ldots \\\left( {H_{N_{t} - 1}\left( {{Bm}_{N_{t} - 1} - {\hat{\overset{\sim}{r}}}_{N_{t} - 1}} \right)} \right)^{T} & & \end{bmatrix}$ $\begin{matrix}{A_{B\;\delta}\overset{\bigtriangleup}{=}\begin{bmatrix}H_{0}^{T} & H_{1}^{T} & \ldots & H_{N_{t} - 1}^{T} & H^{T} & R^{T}\end{bmatrix}} \\{\overset{\bigtriangleup}{=}\begin{bmatrix}H_{T}^{T} & H^{T} & R^{T}\end{bmatrix}}\end{matrix}$$\sum\limits_{H_{T}}^{\;}\;{\overset{\bigtriangleup}{=}{\sigma_{T}^{2}{{blkdiag}\left( {{H_{0}{BH}_{0}^{T}},{H_{1}{BH}_{1}^{T}},\mspace{11mu}\ldots,{H_{N_{t}}{BH}_{N_{t}}^{T}}} \right)}}}$$\sum\limits_{B\;\delta}^{\;}\;{\overset{\bigtriangleup}{=}{\sigma_{T}^{2}{{blkdiag}\left( {\sum\limits_{H_{T}}^{\;}\;{,{\sum\limits_{{\hat{\delta}}^{\bot}}^{\;}\;{,\sum\limits_{R}^{\;}}}}}\; \right)}}}$enables writing the estimator for Bδ as:

$\begin{matrix}{{= {\left( {A_{B\;\delta}^{T}{\sum\limits_{B\;\delta}^{\dagger}\; A_{B\;\delta}}} \right)^{\dagger}A_{B\;\delta}^{T}{\sum\limits_{B\;\delta}^{\dagger}\; m_{B\;\delta}}}}\begin{matrix}{{{cov}{()}}\overset{\bigtriangleup}{=}\left( {A_{B\;\delta}^{T}{\sum\limits_{B\;\delta}^{\dagger}\; A_{B\;\delta}}} \right)^{\dagger}} \\{= {\left\lbrack {{R^{T}{\sum\limits_{R}^{\dagger}\; R}} + {H^{T}{\sum\limits_{\bot}^{\dagger}\; H}} + {H_{T}^{T}{\sum\limits_{H_{T}}^{\dagger}\; H_{T}}}} \right\rbrack^{\dagger}.}}\end{matrix}} & (4.66)\end{matrix}$

A similar formulation can be applied when the range space measurementsare not used as part of the estimation process. In this case, the rangemeasurements are simply omitted from the WLS formulation, and theposition error covariance is modified accordingly. In the case of lowSNR, performing a few iterations of this process provides a solidestimate even though the assumption in Eq. (4.63) does not holdperfectly at first.

At this point the only information that is left unused to improve theestimate of the clock offsets is the terrestrial measurementsR_(t)Bm_(t). These measurements can be expressed as:R _(t)(Bm _(t)−{tilde over ({circumflex over (r)})}_(t))=R _(t) Bδ+R_(t) Be _(t) +R _(t) e _({tilde over (r)}) _(t) ,  (4.67)This equation is solved iteratively, as the uncertainty in the values ofthe clock offsets is dependent on the uncertainty in the terrestrialsource position, which is again dependent on the uncertainty in theclock offsets. Therefore, it is advantageous first to re-estimate theterrestrial source locations using the refined clock offsets. It issimpler to consider the original measurement for the second iteration:Bm _(t)−{tilde over ({circumflex over (r)})}_(t) =Bδ+Be _(t) +e_({tilde over (r)}) _(t) ,  (4.68)It is a good approximation to say that these measurements are relativelyuncorrelated to the clock offset estimations derived from the satellitessignals, Since the little correlation that exists originates from thecurrent position error, which is only partially correlated to the clockoffset measurements derived from the satellites signals.

An expression may be developed as follows for the covariance matrix of:m _(T)

[(Bm ₀−{tilde over ({circumflex over (r)})}₀)^(T), (Bm ₁−{tilde over({circumflex over (r)})}₁)^(T), . . . , (Bm _(N) _(T) ⁻¹−{tilde over({circumflex over (r)})}_(N) _(T) ⁻¹)^(T)]^(T).  (4.69)For this purpose, the error element may be expressed as follows:

$\begin{matrix}\begin{matrix}{{{Be}_{t} + e_{{\overset{\sim}{r}}_{t}}} = {{Be}_{t} + {L_{t}e_{p_{t}}}}} \\{= {{Be}_{t} + {{L_{t}\left( {L_{t}^{T}{\underset{t}{\overset{\dagger}{\sum\limits^{\sim}}}\; L_{t}}} \right)}^{\dagger}L_{t}^{T}{\underset{t}{\overset{\dagger}{\sum\limits^{\sim}}}\;{L_{t}\left( {{Be}_{t} +} \right)}}}}} \\{\overset{\bigtriangleup}{=}{{Be}_{t} + {{\overset{\sim}{L}}_{t}\left( {{Be}_{t} +} \right)}}} \\{= {{\left( {I + {\overset{\sim}{L}}_{t}} \right){Be}_{t}} + {{\overset{\sim}{L}}_{t}.}}}\end{matrix} & (4.70)\end{matrix}$Now we can write the covariance matrix of m_(T):cov(m _(T))

Σ_(T)(Σ_(T))_(i,i)=σ_(T) ²(I+{tilde over (L)} _(i))B(I+{tilde over (L)}_(i))^(T) +{tilde over (L)} _(i)Σ_(Bδ) {tilde over (L)} _(i) ^(T)(Σ_(T))_(i,j) ={tilde over (L)} _(i)Σ_(Bδ){tilde over (L)}_(j)^(T).  (4.71)Combining all measurements of the sensor clock offsets finally gives:

$\begin{bmatrix}m_{T} \\{\hat{\delta}}^{\bot} \\{\hat{\delta}}^{R}\end{bmatrix} = {{\begin{bmatrix}{1_{N_{T}} \otimes I} \\H \\R\end{bmatrix}B\;\delta} + n_{T}}$$\sum\limits_{n_{T}}^{\;}\;{\overset{\bigtriangleup}{=}{{E\left( {n_{T}n_{T}^{T}} \right)} = {{blkdiag}\left( {\sum\limits_{T}^{\;}{,{\sum\limits_{\bot}^{\;}{,\sum\limits_{R}^{\;}}}}} \right)}}}$

The optimal estimator is established using WLS. Defining the following:y _(ƒ)

[m _(T) ^(T), {circumflex over (δ)}^(⊥) ^(T) , {circumflex over (δ)}^(R)^(T) ]^(T)A _(ƒ)

[1_(N) _(T)

I)^(T) , H ^(T) , R ^(T)]^(T)Σ_(ƒ)

blkdiag(Σ_(T), Σ_(⊥), Σ_(R))the estimator is given by:

=A _(ƒ) ^(T)Σ_(ƒ) ^(†) A _(ƒ))^(†) A _(ƒ) ^(T)Σ_(ƒ) ^(†) y _(ƒ)cov(

)=(A _(ƒ) ^(T)Σ_(ƒ) ^(†) A _(ƒ))^(†).

After estimating the clock offsets, controller 28 may again re-estimatethe locations of terrestrial sources 24, and repeat until convergence asdescribed above.

In some embodiments it is possible to skip the steps that use RBm_(t)and return, or alternately do the final refinement using gradientmethods on the original cost function and measurement vector.

Based on the above estimators, it can be shown that even with only onesatellite source 26, it is possible to achieve accurate synchronizationand localization as long as there is a terrestrial source 24 that isreceived by at least six or seven sensors (depending on whether they aredistributed in a plane or in 3D space). In other words, for anycombination of satellite sources 26 and terrestrial sources 24, as longas they are received by sufficient number of sensors 22, the estimatorsdefined above are well posed and have performance directly related tothe measurement error.

The Navigation Problem

The problem of navigation of a moving object in system 20 is equivalentto assuming that a sensor of unknown position is added to the network.This addition does not change the form of the measurement model.Considering only celestial sources, for example, the only difference isthat now one row of the G matrix is unknown, and another clock offset isadded.

Defining s^(T) as the unknown sensor position and δ_(s) as its clockoffset, the estimator may be written as:

$\begin{matrix}\begin{matrix}{{Bm}_{t} = {{{- {B\begin{bmatrix}G \\s^{T}\end{bmatrix}}}{\overset{\_}{p}}_{t}} + {B\begin{bmatrix}\delta \\\delta_{S}\end{bmatrix}} + {Be}_{t}}} \\{= {{{- {B\left( {\begin{bmatrix}G \\0\end{bmatrix} + \begin{bmatrix}0 \\s^{T}\end{bmatrix}} \right)}}{\overset{\_}{p}}_{t}} + {B\begin{bmatrix}\delta \\0\end{bmatrix}} + {B\begin{bmatrix}0 \\\delta_{S}\end{bmatrix}} + {Be}_{t}}}\end{matrix} & (4.72)\end{matrix}$Assuming 1^(T)δ=0, as is expected from the estimation process of theclock offsets, leads to:

${Bm}_{t} = {{{- {B\begin{bmatrix}G \\0\end{bmatrix}}}{\overset{\_}{p}}_{t}} + {{B\begin{bmatrix}0 \\s^{T}\end{bmatrix}}{\overset{\_}{p}}_{t}} + \begin{bmatrix}\delta \\0\end{bmatrix} + {B\begin{bmatrix}0 \\\delta_{S}\end{bmatrix}} + {Be}_{t}}$Rearranging the equation gives:

${{Bm}_{t} = {{\Gamma_{t}\left( {{\overset{\_}{p}}_{t},\delta} \right)} + {{A_{t}\left( {\overset{\_}{p}}_{t} \right)}\begin{bmatrix}s \\\delta_{S}\end{bmatrix}} + {Be}_{t}}},$wherein:

$\begin{matrix}{{{\Gamma_{t}\left( {{\overset{\_}{p}}_{t},\delta} \right)}\overset{\bigtriangleup}{=}{{{- {B\begin{bmatrix}G \\0\end{bmatrix}}}{\overset{\_}{p}}_{t}} + \begin{bmatrix}\delta \\0\end{bmatrix}}}{{A_{t}\left( {\overset{\_}{p}}_{t} \right)}\overset{\bigtriangleup}{=}{- {{B\begin{bmatrix}0 & 0 \\{\overset{\_}{p}}_{t}^{T} & {- 1}\end{bmatrix}}.}}}} & (4.73)\end{matrix}$

Collecting the equations for all t gives:

${m = {{\Gamma\left( {\left\{ {\overset{\_}{p}}_{i} \right\}_{i = 0}^{N_{S} - 1},\delta} \right)} + {{A\left( \left\{ {\overset{\_}{p}}_{i} \right\}_{i = 0}^{N_{S} - 1} \right)}\begin{bmatrix}s \\\delta_{S}\end{bmatrix}} + e}},$wherein:Γ({ p _(i)}_(i=0) ^(N) ^(s) ⁻¹, δ)

[Γ₀ ^(T), Γ₁ ^(T), . . . , Γ_(N) _(s) ⁻¹]^(T)A({ p _(i)}_(i=0) ^(N) ^(s) ⁻¹)

[A ₀ ^(T) , A ₁ ^(T) , . . . , A _(N) _(s) ⁻¹]^(T).  (4.74)

Substituting for the satellite directions and the sensor clock offsetswith their estimates (based on the methods described above without theuse of the sensor of unknown position) gives:

${m = {{\Gamma\left( {\left\{ {\hat{\overset{\_}{p}}}_{i} \right\}_{i = 0}^{N_{S} - 1},\hat{\delta}} \right)} + {{A\left( \left\{ {\hat{\overset{\_}{p}}}_{i} \right\}_{i = 0}^{N_{S} - 1} \right)}\begin{bmatrix}s \\\delta_{S}\end{bmatrix}} + e}},$The navigation problem is thus reduced to a weighed least squaresestimation of the values:

$\begin{matrix}{\begin{bmatrix}\hat{s} \\{\hat{\delta}}_{S}\end{bmatrix} = {\left( {A^{T}A} \right)^{- 1}{{A^{T}\left( {m - \Gamma} \right)}.}}} & (4.75)\end{matrix}$

This estimation can be used as an initial point for the minimization ofthe full ML cost presented in Eq. (4.13), with the only difference beingthat the position of one of the sensors is unknown. Starting from theinitial estimate given by Eq. (4.75), controller 28 can use gradientmethods to find the optimal point in terms of ML cost.

The problem framed above can be extended to the case of multiple sensorsof unknown position by solving the LS estimation for each of thesesensors. The ML estimator can use these initial estimates to find theoptimum of the full ML cost, and in so doing also refine the estimatesof the known sensors clock offsets and the satellite directions.Terrestrial sources can also be incorporated into the ML cost function.

APPLICATION EXAMPLES

FIG. 5 is a schematic, pictorial illustration showing a navigationsystem 70 based on a network of sensors 22, in accordance with anembodiment of the present invention. The sensor network in this exampleis deployed in an airport, with controller 28 housed in a control center72. Airplanes 74 are also allowed to carry sensors 76 (or to be able totransmit or reflect a signal that sensors 22 can detect). Sources thatcan be used by the system include ground sources 24 (such as beacons,radar signals), aerial sources (airplane beacons 76), and sky sources(satellites 26). System 70 will be generally immune to jamming andspoofing, unlike GNSS systems.

There are two possible modes of navigation in system 70:

a. Standalone mode: The measurements of sensor 76 are used only toresolve its location and time offset. This mode is linear in the sensorposition and clock offset, and therefore fast, and could be carried outby a processor (not shown) on board airplane 74.

b. Joint mode: The measurements of sensor 76 are also used to improveall other system estimates, including the positions and time estimatesof other airplane sensor.

The arrangement shown in FIG. 5 enables accurate airplane navigation,and the measurements of airplane sensors 76 can also be used to improvesystem accuracy and integrity even further using the methods describedabove. As noted earlier, the readings of sensor 76 can be used toimprove clock synchronization, terrestrial source localization andsatellite direction accuracy. Navigation is still possible even if onlyterrestrial sources are present, although the processing required istypically more complex.

FIG. 6 is a schematic, pictorial illustration of a computer system 80with sensor-based synchronization, in accordance with an embodiment ofthe present invention. System 80 may comprise, for example, acomputerized securities trading system, in which servers 82 areconnected by a network 84 and user sensors 22 in accuratelysynchronizing their respective clocks based on signals from satellitesources 26.

Trading time is of great importance in such systems, and time-stampingis regulated by law. In systems that rely on GPS-based synchronization,if the GPS receiver is jammed and accurate time stamping is notavailable, trading must stop until it is recovered. Spoofing of thetrading system time could cause trades to be recorded at a loss.

In contrast to existing solutions, servers 82 use sensors 22 toimplement spoof-proof and jamming-proof time-stamping, based on themethods of clock synchronization that are described above. As timetransfer is enabled without the use of GNSS, the absolute time can bemaintained by all servers 82, which are unlikely to be GNSS jammed atthe same time, or the servers as a group could perform time holdoveruntil the total jamming has passed. System 80 is thus able to maintainnormal operation even under jamming attacks and to detect and rejectattempts at source spoofing quickly and robustly.

FIG. 7 is a schematic, pictorial illustration of a communication system90 with sensor-based synchronization, in accordance with an embodimentof the present invention. System 90 comprises multiple base stations(BS) 92, which are coupled via a network 94 to controller 28. The basestations collect TOA information gathered by user equipment 96, andcontroller 28 (or distributed controllers in base stations 92) uses thisinformation in synchronizing the base station clocks. System 90 differsfrom the embodiments described above in that the positions of thesources (base stations 92) are known precisely, while the positions ofthe sensors (user equipment 96) may not be known accurately. Clocksynchronization is accomplished, however, substantially as describedabove.

The precise synchronization that is achieved in system 90 can be used invarious ways to improve network performance, such as reducinginterference, increasing network capacity or bits per Hz, improvingresource allocation, enabling Cooperative Multipoint (COMP) transmissionin small and macrocell deployments, and supporting time-basedlocalization techniques such as OTDOA (observed time difference ofarrival), and UTDOA.

Assuming system 90 operates in accordance with the Long-Term Evolution(LTE) family of standards, clock synchronization in the system can beimplemented as follows:

-   -   1. Collect OTDOA measurements or equivalent type of measurements        from user equipment 96.        -   a. Such measurements may be invoked by enabling active            location request generation, which generates localization            requests to the user equipment in the area of the base            stations that are to be synchronized.        -   b. If the user equipment also calculates its GPS position,            the base stations may retrieve this information, as well,            since then the UE will be regarded as a source (as well as a            sensor) of known position and can improve performance as            described above.    -   2. Process all measurements in controller 28, as described        above, and return the user equipment positions and base stations        clock offsets and, if needed, their frequency offsets.    -   3. Update the base stations clocks accordingly.        The rate at which the process is repeated is dependent on the        stability of the base stations clocks relative to one another.        Optionally, the rate at which the above process is performed can        be reduced by modeling the dynamics of the changing clock        offsets.

Although the above example, relates specifically to LTE, other mobilenetwork technologies (such as UMTS) that provide OTDOA-like measurementscan use the methods described herein.

One such measurement type is UTDOA (Uplink TDOA), which is performed bydedicated sensors called LMUs (Location Measurement Units). The LMU canbe an integral part of the base station or a separate unit, installed onthe base station or in another position. When using UTDOA measurements,the emitters are the user equipment (UE), and the LMUs are consideredthe sensors of known position. Using the methods described herein, theLMUs clocks are synchronized, while the UE locations are also resolved.Moreover, if the LMU clock is correlated to the base station clock, thenthe base stations and LMUs are also synchronized, enabling all of thebenefits stated above and eliminating the need for a dedicated externalsynchronization unit in each LMU.

In another embodiment (not shown explicitly in the figures), sensors areused to gather information on the environment, and the locations of thesensors need to be determined (as in the case of Real Time LocationServices, for example). These environmental sensors may be regarded as“sources” in the context of this embodiment, and localization of thesesensors may be achieved based solely on the signals that they emit. Thesensor measurements of these signals constitute the information gatheredby the sensors that is relayed to a central node, so that no specialtransmissions are required for achieving sensor localization.

In such a network, in other words, there are two types of sensors:

-   1. Localization sensors 22 (which measure TOA), preferably of known    position.-   2. General sensing sensors, which gather sensing information on    their environment (such as temperature, vibration, humidity, etc.)    and serve as sources 24.    The synchronization that is achieved among the type 1 sensors can be    used to efficiently manage the network.

In some embodiments, synchronization may be the major interest and notlocalization. This situation can arise when, for example, a wirelessnetwork needs to manage wakeup times of its nodes. An embodiment of thepresent invention provides a good solution for such cases, sinceemitters are usually present in any environment. Thus, even if onlypairwise distances of the sensors are known, synchronization can beachieved to some extent.

Additionally or alternatively, the techniques described herein can beused to check whether an existing system is synchronized and calibratedproperly with respect to time offsets and skews, i.e., to detect theoccurrence of synchronization faults. If the clock offsets and skewsderived using the methods described above do not match those in theexisting system, there may be a malfunction in the hardware or softwareof the existing system. In this case, the time offset and skew of thesystem can be recalibrated based on the synchronization described by theabove methods, and this synchronization can also be used to provide abackup clock or even a set of such clocks for the system, thus affordingmaximal availability and reliability.

Although the embodiments disclosed above relate specifically to problemsof sensor synchronization and emitter localization (along withnavigation and spoofing detection), other problems that can beformulated in an equivalent manner can be also solved by the methodsthat are described herein. For example, the principles of the techniquesdescribed in the present patent application may be extended toequivalent embodiments that arise from the time/frequency dualitybetween DTOA measurements and differential Doppler measurements. Allsuch solutions are considered to be within the scope of the presentinvention.

It will thus be appreciated that the embodiments described above arecited by way of example, and that the present invention is not limitedto what has been particularly shown and described hereinabove. Rather,the scope of the present invention includes both combinations andsubcombinations of the various features described hereinabove, as wellas variations and modifications thereof which would occur to personsskilled in the art upon reading the foregoing description and which arenot disclosed in the prior art.

The invention claimed is:
 1. A method for transmitter operation,comprising: receiving recorded time of arrival (TOA) measurements oftransmitted signals from a plurality of sensors, at unknown locations,the signals being transmitted from at least a group of sources in anetwork of sources, which comprise respective clocks that are notmutually synchronized, and which are configured to transmit the signalsat times determined according to the respective clocks, wherein thesensors do not serve as sources of the transmitted signals; obtaininglocation information comprising respective source locations of thesources; calculating synchronization values for the sources comprisingat least one of respective offsets and skews of the clocks thereof as afunction of the location information comprising the respective sourcelocations, and the received recorded TOA measurements from the sensors,without dependence on measurements of transmissions emitted by thesensors to the sources; and synchronizing the respective clocks of thesources or finding locations of the sensors based on the calculatedsynchronization values; wherein calculating the synchronization valuescomprises applying an estimator to a set of equations relating therecorded TOA measurements, source synchronization values and the sourceand sensor locations; wherein applying the estimator comprises applyingan iterative optimization process to the set of the equations; whereinthe iterative optimization process iteratively calculate estimatedlocations of the plurality of sensors and intermediate synchronizationvalues while taking into account a random measurement error and therespective source locations of the sources.
 2. The method according toclaim 1, wherein the network of sources comprises base stations in acellular communications network, and the sensors comprise user equipmentin the cellular communications network, and wherein calculating thesynchronization values comprises calculating the synchronization valuesas functions of signals received by a plurality of the user equipments.3. A network system, comprising: a network of sources, which compriserespective clocks that are not mutually synchronized, and which areconfigured to transmit, at times determined according to the respectiveclocks, from at least a group of the sources, respective signals; and aprocessor, which is configured to receive recorded time of arrival (TOA)measurements of the transmitted signals from a plurality of sensors, atunknown locations, to process the recorded TOA measurements, usinglocation information comprising respective source locations of thesources, so as to calculate synchronization values for the sourcescomprising at least one of respective offsets and skews of the clocksthereof as a function of the location information comprising therespective source locations, and the received recorded TOA measurementsfrom the sensors, without dependence on measurements of transmissionsemitted by the sensors to the sources, and to synchronize the respectiveclocks of the sources based on the calculated synchronization values;wherein the synchronization values are calculated by applying anestimator to a set of equations relating the recorded TOA measurements,source synchronization values and the source and sensor locations;wherein applying the estimator comprises applying an iterativeoptimization process to the set of the equations; wherein the iterativeoptimization process iteratively calculate estimated locations of theplurality of sensors and intermediate synchronization values whiletaking into account a random measurement error and the respective sourcelocations of the sources.
 4. The system according to claim 3, whereinthe network of sources comprises base stations in a cellularcommunications network, and the sensors comprise user equipment in thecellular communications network, and wherein the processor is configuredto calculate synchronization values as functions of signals received bya plurality of the user equipment.
 5. The system according to claim 4,wherein the processor is configured to resolve respective equipmentlocations of the user equipment based on the signals received by theuser equipment.
 6. A a non-transitory computer-readable medium in whichprogram instructions are stored, which instructions, when read by acomputer, cause the computer to receive respective time of arrival (TOA)measurements of signals received by a plurality of sensors, at unknownlocations, from at least a group of sources in a network of sources,which comprise respective clocks that are not mutually synchronized, andwhich are configured to transmit the signals at times determinedaccording to the respective clocks, and to process the recorded TOAmeasurements, using location information comprising respective sourcelocations of the sources, so as to calculate synchronization values forthe sources comprising at least one of respective offsets and skews ofthe clocks thereof as a function of the location information comprisingthe respective source locations and the received recorded TOAmeasurements from the sensors, without dependence on measurements oftransmissions emitted by the sensors to the sources, and to synchronizethe respective clocks of the sources based on the calculatedsynchronization values; wherein the synchronization values arecalculated by applying an estimator to a set of equations relating therecorded TOA measurements, source synchronization values and the sourceand sensor locations; wherein applying the estimator comprises applyingan iterative optimization process to the set of the equations; whereinthe iterative optimization process iteratively calculate estimatedlocations of the plurality of sensors and intermediate synchronizationvalues while taking into account a random measurement error and therespective source locations of the sources.
 7. The method according toclaim 1, wherein synchronizing the respective clocks comprisesestimating offsets and skews between the respective clocks.
 8. Themethod according to claim 1, wherein the method comprises computing thesensor locations based on the source locations and the recorded TOAmeasurements.
 9. The method according to claim 1, wherein applying theiterative optimization process comprises performing a convexoptimization.
 10. The method according to claim 1, and comprisingdetecting a fault in the network based on the location information andthe recorded TOA measurements.
 11. The system according to claim 4,wherein the cellular communications network operates in accordance witha Long-Term Evolution (LTE) standard family, and wherein the userequipment is configured to make observed time difference of arrival(OTDOA) measurements based on the signals received from the basestations, and the processor is configured to perform at least one ofsynchronizing the operation of the base stations and finding locationsof the user equipment based on the OTDOA measurements.
 12. The systemaccording to claim 3, wherein the processor is configured to estimateoffsets and skews between the respective clocks.
 13. The systemaccording to claim 3, wherein the processor is configured to compute thesensor locations based on the source locations and the recorded TOAmeasurements.
 14. The system according to claim 13, wherein theprocessor is configured to apply an estimator to a set of equationsrelating the recorded TOA measurements and the source and sensorlocations in order to compute the sensor locations.
 15. The systemaccording to claim 14, wherein the processor is configured to apply theestimator to the set of the equations in an iterative optimizationprocess.
 16. The system according to claim 15, wherein the iterativeoptimization process comprises a convex optimization.
 17. The systemaccording to claim 3, wherein the processor is configured to detect afault in the network based on the location information and the recordedTOA measurements.
 18. The system according to claim 6, wherein theinstructions cause the computer to compute the sensor locations based onthe source locations and the recorded TOA measurements.
 19. The methodaccording to claim 2, comprising transmitting a localization request tosensors in an area of the sources to cause the user equipment to providethe recorded TOA measurements.
 20. The system according to claim 3,wherein the TOA measurements comprise differential time of arrivalmeasurements.
 21. The system according to claim 3, wherein the TOAmeasurements comprise times of arrival.
 22. The method according toclaim 2, comprising resolving respective equipment locations of the userequipment based on the calculated synchronization values.